From 2610260f56f43c5f6531742890e4a55221f22222 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 10 Feb 2011 18:56:37 +0000 Subject: [PATCH] Double-indexing for energy grid, made Lists and Dictionaries more generic with interfaces. Sampling of nuclides. --- ChangeLog | 31 ++ src/Makefile | 11 +- src/ace.f90 | 89 ++-- src/cross_section.f90 | 66 ++- src/data_structures.f90 | 985 +++++++++++++++++++++++++++++++--------- src/energy_grid.f90 | 86 ++-- src/fileio.f90 | 275 ++++++++--- src/geometry.f90 | 323 ++++++++----- src/global.f90 | 29 +- src/input_sample | 14 +- src/main.f90 | 41 +- src/physics.f90 | 161 ++++++- src/search.f90 | 16 +- src/source.f90 | 4 +- src/string.f90 | 6 +- src/types.f90 | 177 +++++--- 16 files changed, 1711 insertions(+), 603 deletions(-) diff --git a/ChangeLog b/ChangeLog index 80e5465e1a..5cbfb09c18 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,34 @@ +2011-02-09 Paul Romano + + * ace.f90: Fixed bug in the way xs_continuous was allocated. If, + e.g., 92235.03c was listed twice in a file, two spots would be + allocated in the array instead of one. + * cross_section.f90: Added calculation of total material + cross-section. + * data_structures.f90: Added list routines for each different + type (integer, real, key/value pairs) and an interface for + them. Dictionaries can now take character keys or integer + keys. New hash function for integers is simple mod function. + * energy_grid.f90: Removed grid_count since we now have a generic + list_size routine in data_structures. Added original_indices + routine in order to do SERPENT-style double-indexing. + * fileio.f90: Fixed normalize_ao routine to include calculations + based on atom or gram density (this should probably be moved to + different module). Added get_next_line routine that reads a full + line of input including any continuation lines marked by an + ampersand. + * global.f90: Added separate dictionary types (CI and II), added + pointers for current cell, surface, and material. Added e_grid and + n_grid + * physics.f90: Neutron energy sampled at beginning of + transport. Added watt_spectrum and maxwell_spectrum + routines. Changed collision routine so that actual nuclide is + sampled. + * search.f90: Fixed binary_search routine + * string.f90: Changed concatenate subroutine into a function + * types.f90: Moved List and Dictionary types from data_structures + to here (including new types). LinkedListGrid is now ListReal. + 2011-02-03 Paul Romano * ace.f90: Added initial routines for reading continuous-energy diff --git a/src/Makefile b/src/Makefile index bb1eb4218f..341826cf18 100644 --- a/src/Makefile +++ b/src/Makefile @@ -15,7 +15,7 @@ modules = ace.f90 \ search.f90 \ source.f90 \ string.f90 \ - types.f90 \ + types.f90 main_objects = $(modules:.f90=.o) $(main:.f90=.o) test_objects = $(modules:.f90=.o) $(test:.f90=.o) @@ -50,16 +50,17 @@ unittest: $(test_objects) ace.o: global.o output.o string.o fileio.o string.o cross_section.o: global.o string.o data_structures.o output.o data_structures.o: global.o -energy_grid.o: global.o +energy_grid.o: global.o output.o data_structures.o fileio.o: types.o global.o string.o output.o data_structures.o -geometry.o: types.o global.o output.o string.o +geometry.o: types.o global.o output.o string.o data_structures.o global.o: types.o main.o: global.o fileio.o output.o geometry.o mcnp_random.o \ source.o physics.o cross_section.o data_structures.o \ ace.o energy_grid.o output.o: global.o -physics.o: types.o global.o mcnp_random.o geometry.o output.o -search.o: +physics.o: types.o global.o mcnp_random.o geometry.o output.o \ + search.o +search.o: output.o source.o: global.o mcnp_random.o string.o: global.o output.o types.o: diff --git a/src/ace.f90 b/src/ace.f90 index 0227ce2216..7a5bce9e32 100644 --- a/src/ace.f90 +++ b/src/ace.f90 @@ -32,69 +32,60 @@ contains type(AceThermal), pointer :: ace_thermal => null() integer :: i, j integer :: index - character(10) :: table + character(10) :: key character(250) :: msg integer :: n - integer :: n_continuous - integer :: n_thermal integer :: index_continuous integer :: index_thermal + type(ListKeyValueCI), pointer :: elem + + call dict_create(ace_dict) ! determine how many continuous-energy tables and how many S(a,b) ! thermal scattering tables there are - n_continuous = 0 - n_thermal = 0 - do i = 1, n_materials - mat => materials(i) - do j = 1, mat%n_isotopes - index = mat%isotopes(j) - table = xsdatas(index)%id - n = len_trim(table) - call lower_case(table) - select case (table(n:n)) - case ('c') - n_continuous = n_continuous + 1 - case ('t') - n_thermal = n_thermal + 1 - case default - msg = "Unknown cross section table type: " // table - call error(msg) - end select - end do - end do - - ! allocate arrays for ACE table storage - allocate(xs_continuous(n_continuous)) - allocate(xs_thermal(n_thermal)) - - ! loop over all materials - call dict_create(ace_dict) index_continuous = 0 index_thermal = 0 do i = 1, n_materials mat => materials(i) do j = 1, mat%n_isotopes index = mat%isotopes(j) - table = xsdatas(index)%id - n = len_trim(table) - call lower_case(table) - select case (table(n:n)) + key = xsdatas(index)%id + n = len_trim(key) + call lower_case(key) + select case (key(n:n)) case ('c') - if (.not. dict_has_key(ace_dict, table)) then + if (.not. dict_has_key(ace_dict, key)) then index_continuous = index_continuous + 1 - call dict_add_key(ace_dict, table, index_continuous) - call read_ACE_continuous(index_continuous, index) + call dict_add_key(ace_dict, key, index_continuous) mat%table(j) = index_continuous + else + mat%table(j) = dict_get_key(ace_dict, key) end if case ('t') - index_thermal = index_thermal + 1 n_thermal = n_thermal + 1 case default - msg = "Unknown cross section table type: " // table + msg = "Unknown cross section table type: " // key call error(msg) end select end do end do + + n_continuous = index_continuous + n_thermal = index_thermal + + ! allocate arrays for ACE table storage + allocate(xs_continuous(n_continuous)) + allocate(xs_thermal(n_thermal)) + + ! loop over all nuclides in xsdata + index_continuous = 0 + do i = 1, size(xsdatas) + key = xsdatas(i)%alias + if (dict_has_key(ace_dict, key)) then + index_continuous = index_continuous + 1 + call read_ACE_continuous(index_continuous, i) + end if + end do end subroutine read_xs @@ -130,14 +121,14 @@ contains table => xs_continuous(index_table) ! Check if input file exists and is readable - inquire( FILE=filename, EXIST=file_exists, READ=readable ) - if ( .not. file_exists ) then + inquire(FILE=filename, EXIST=file_exists, READ=readable) + if (.not. file_exists) then msg = "ACE library '" // trim(filename) // "' does not exist!" - call error( msg ) - elseif ( readable(1:3) /= 'YES' ) then + call error(msg) + elseif (readable(1:3) == 'NO') then msg = "ACE library '" // trim(filename) // "' is not readable! & &Change file permissions with chmod command." - call error( msg ) + call error(msg) end if ! display message @@ -145,7 +136,12 @@ contains call message(msg, 6) ! open file - open(file=filename, unit=in, status='old', action='read') + open(file=filename, unit=in, status='old', & + & action='read', iostat=ioError) + if (ioError /= 0) then + msg = "Error while opening file: " // filename + call error(msg) + end if found_xs = .false. do while (.not. found_xs) @@ -217,6 +213,7 @@ contains ! determine number of energy points NE = NXS(3) + table%n_grid = NE ! allocate storage for arrays allocate(table%energy(NE)) @@ -319,7 +316,7 @@ contains ! cross section table !===================================================================== - subroutine read_ACE_thermal( filename ) + subroutine read_ACE_thermal(filename) character(*), intent(in) :: filename diff --git a/src/cross_section.f90 b/src/cross_section.f90 index 79d629d938..4eafe7cafc 100644 --- a/src/cross_section.f90 +++ b/src/cross_section.f90 @@ -1,9 +1,9 @@ module cross_section + use global use string, only: split_string, str_to_real - use global, only: str_to_int, max_words, xsdata_dict, xsdatas - use data_structures, only: Dictionary, dict_create, dict_add_key, dict_get_key - use output, only: error + use data_structures, only: dict_create, dict_add_key, dict_get_key + use output, only: error, message use types, only: xsData implicit none @@ -31,6 +31,9 @@ contains integer :: index integer :: ioError + msg = "Reading cross-section summary file..." + call message(msg, 5) + ! Construct filename n = len(path) if (path(n:n) == '/') then @@ -44,14 +47,19 @@ contains if ( .not. file_exists ) then msg = "Cross section summary '" // trim(filename) // "' does not exist!" call error( msg ) - elseif ( readable(1:3) /= 'YES' ) then + elseif ( readable(1:3) == 'NO' ) then msg = "Cross section summary '" // trim(filename) // "' is not readable!" & & // "Change file permissions with chmod command." call error( msg ) end if - ! open file - open(file=filename, unit=in, status='old', action='read') + ! open xsdata file + open(file=filename, unit=in, status='old', & + & action='read', iostat=ioError) + if (ioError /= 0) then + msg = "Error while opening file: " // filename + call error(msg) + end if ! determine how many lines count = 0 @@ -110,4 +118,50 @@ contains end subroutine read_xsdata +!===================================================================== +! MATERIAL_TOTAL_XS calculates the total macroscopic cross section of +! each material for use in sampling path lengths +!===================================================================== + + subroutine material_total_xs() + + type(Material), pointer :: mat => null() + type(AceContinuous), pointer :: acecont => null() + integer :: i, j, k + integer :: index + real(8) :: xs + real(8) :: density + character(250) :: msg + + msg = "Creating material total cross-sections..." + call message(msg, 4) + + ! loop over all materials + do i = 1, n_materials + + ! allocate storage for total xs + mat => materials(i) + density = mat%atom_density + allocate(mat%total_xs(n_grid)) + + ! loop over each energy on grid + do j = 1, n_grid + + mat%total_xs(j) = 0.0_8 + + ! loop over each isotope in material + do k = 1, mat%n_isotopes + ! find ACE cross section table for this isotope + index = mat%table(k) + acecont => xs_continuous(index) + + ! add contribution to total cross-section + mat%total_xs(j) = mat%total_xs(j) + density * & + mat%atom_percent(j) * acecont%sigma_t(index) + end do + end do + end do + + end subroutine material_total_xs + end module cross_section diff --git a/src/data_structures.f90 b/src/data_structures.f90 index 1882f570b0..7462d793bb 100644 --- a/src/data_structures.f90 +++ b/src/data_structures.f90 @@ -14,138 +14,321 @@ module data_structures ! to ListData, dict_create, dict_add_key, and dict_get_key). !===================================================================== - use types, only: ListData, LinkedList, HashList, Dictionary + use types implicit none - ! Hide objects that do not need to be publicly accessible - private :: ListData - private :: HashList - private :: LinkedList - private :: list_create - private :: list_destroy - private :: list_count - private :: list_next - private :: list_insert - private :: list_insert_head - private :: list_delete_element - private :: list_get_data - private :: list_put_data - private :: dict_get_elem - private :: dict_hashkey - integer, parameter, private :: hash_size = 4993 integer, parameter, private :: multiplier = 31 integer, parameter :: DICT_NULL = 0 +!===================================================================== +! LIST Interfaces -- these allow one to use a single subroutine or +! function for various types of lists, namely those consisting of +! integers, reals, or (key,value) pairs. +!===================================================================== + + interface list_create + module procedure list_real_create, list_int_create + module procedure list_kvci_create, list_kvii_create + end interface + interface list_insert + module procedure list_int_insert, list_real_insert + module procedure list_kvci_insert, list_kvii_insert + end interface + interface list_delete + module procedure list_int_delete, list_real_delete + module procedure list_kvci_delete, list_kvii_delete + end interface + interface list_size + module procedure list_int_size, list_real_size + module procedure list_kvci_size, list_kvii_size + end interface + interface list_insert_head + module procedure list_int_insert_head, list_real_insert_head + module procedure list_kvci_insert_head, list_kvii_insert_head + end interface + interface list_delete_element + module procedure list_int_del_element, list_real_del_element + module procedure list_kvci_del_element, list_kvii_del_element + end interface + +!===================================================================== +! DICTIONARY Interfaces -- these allow one to use a single subroutine +! or function for dictionaries with varying types of (key,value) pairs +!===================================================================== + + interface dict_create + module procedure dict_ci_create, dict_ii_create + end interface + interface dict_delete + module procedure dict_ci_delete, dict_ii_delete + end interface + interface dict_add_key + module procedure dict_ci_add_key, dict_ii_add_key + end interface + interface dict_delete_key + module procedure dict_ci_delete_key, dict_ii_delete_key + end interface + interface dict_get_key + module procedure dict_ci_get_key, dict_ii_get_key + end interface + interface dict_has_key + module procedure dict_ci_has_key, dict_ii_has_key + end interface + interface dict_get_elem + module procedure dict_ci_get_elem, dict_ii_get_elem + end interface + interface dict_hashkey + module procedure dict_ci_hashkey, dict_ii_hashkey + end interface + contains !===================================================================== -! LIST_CREATE creates and initializes a list -! Arguments: -! list Pointer to new linked list -! data The data for the first element -! Note: -! This version assumes a shallow copy is enough -! (that is, there are no pointers within the data -! to be stored) -! It also assumes the argument list does not already -! refer to a list. Use list_destroy first to -! destroy up an old list. +! LIST_INT_CREATE creates and initializes a list of integers !===================================================================== - subroutine list_create( list, data ) + subroutine list_int_create(list, data) - type(LinkedList), pointer :: list - type(ListData), intent(in) :: data + type(ListInt), pointer :: list + integer, intent(in) :: data - allocate( list ) + allocate(list) list%next => null() list%data = data - end subroutine list_create + end subroutine list_int_create !===================================================================== -! LIST_DESTROY destroys an entire list -! Arguments: -! list Pointer to the list to be destroyed -! Note: -! This version assumes that there are no -! pointers within the data that need deallocation +! LIST_REAL_CREATE creates and initializes a list of real(8)s !===================================================================== - subroutine list_destroy( list ) + subroutine list_real_create(list, data) - type(LinkedList), pointer :: list + type(ListReal), pointer :: list + real(8), intent(in) :: data - type(LinkedList), pointer :: current - type(LinkedList), pointer :: next + allocate(list) + list%next => null() + list%data = data + + end subroutine list_real_create + +!===================================================================== +! LIST_KVCI_CREATE creates and initializes a list of (key,value) pairs +!===================================================================== + + subroutine list_kvci_create(list, data) + + type(ListKeyValueCI), pointer :: list + type(KeyValueCI), intent(in) :: data + + allocate(list) + list%next => null() + list%data = data + + end subroutine list_kvci_create + +!===================================================================== +! LIST_KVII_CREATE creates and initializes a list of (key,value) pairs +!===================================================================== + + subroutine list_kvii_create(list, data) + + type(ListKeyValueII), pointer :: list + type(KeyValueII), intent(in) :: data + + allocate(list) + list%next => null() + list%data = data + + end subroutine list_kvii_create + +!===================================================================== +! LIST_INT_DELETE deallocates all data in a list of integers +!===================================================================== + + subroutine list_int_delete(list) + + type(ListInt), pointer :: list + type(ListInt), pointer :: current + type(ListInt), pointer :: next current => list - do while ( associated(current%next) ) + do while (associated(current%next)) next => current%next - deallocate( current ) + deallocate(current) current => next enddo - end subroutine list_destroy + end subroutine list_int_delete !===================================================================== -! LIST_COUNT count the number of items in the list -! Arguments: -! list Pointer to the list +! LIST_REAL_DELETE deallocates all data in a list of real(8)s !===================================================================== - integer function list_count( list ) + subroutine list_real_delete(list) - type(LinkedList), pointer :: list + type(ListReal), pointer :: list + type(ListReal), pointer :: current + type(ListReal), pointer :: next - type(LinkedList), pointer :: current - type(LinkedList), pointer :: next + current => list + do while (associated(current%next)) + next => current%next + deallocate(current) + current => next + enddo - if ( associated(list) ) then - list_count = 1 + end subroutine list_real_delete + +!===================================================================== +! LIST_KVCI_DELETE deallocates all data in a list of (key,value) pairs +!===================================================================== + + subroutine list_kvci_delete(list) + + type(ListKeyValueCI), pointer :: list + type(ListKeyValueCI), pointer :: current + type(ListKeyValueCI), pointer :: next + + current => list + do while (associated(current%next)) + next => current%next + deallocate(current) + current => next + enddo + + end subroutine list_kvci_delete + +!===================================================================== +! LIST_KVII_DELETE deallocates all data in a list of (key,value) pairs +!===================================================================== + + subroutine list_kvii_delete(list) + + type(ListKeyValueII), pointer :: list + type(ListKeyValueII), pointer :: current + type(ListKeyValueII), pointer :: next + + current => list + do while (associated(current%next)) + next => current%next + deallocate(current) + current => next + enddo + + end subroutine list_kvii_delete + +!===================================================================== +! LIST_INT_SIZE counts the number of items in a list of integers +!===================================================================== + + function list_int_size(list) result(n) + + type(ListInt), pointer :: list + type(ListInt), pointer :: current + integer :: n + + if (associated(list)) then + n = 1 current => list - do while ( associated(current%next) ) + do while (associated(current%next)) current => current%next - list_count = list_count + 1 + n = n + 1 enddo else - list_count = 0 + n = 0 endif - end function list_count + end function list_int_size !===================================================================== -! LIST_NEXT returns the next element (if any) -! Arguments: -! elem Element in the linked list -! Result: +! LIST_REAL_SIZE counts the number of items in a list of real(8)s !===================================================================== - function list_next( elem ) result(next) + function list_real_size(list) result(n) - type(LinkedList), pointer :: elem - type(LinkedList), pointer :: next + type(ListReal), pointer :: list + type(ListReal), pointer :: current + integer :: n - next => elem%next + if (associated(list)) then + n = 1 + current => list + do while (associated(current%next)) + current => current%next + n = n + 1 + enddo + else + n = 0 + endif - end function list_next + end function list_real_size !===================================================================== -! LIST_INSERT inserts a new element +! LIST_KVCI_SIZE counts the number of items in a list of (key,value) +! pairs +!===================================================================== + + function list_kvci_size(list) result(n) + + type(ListKeyValueCI), pointer :: list + type(ListKeyValueCI), pointer :: current + integer :: n + + if (associated(list)) then + n = 1 + current => list + do while (associated(current%next)) + current => current%next + n = n + 1 + enddo + else + n = 0 + endif + + end function list_kvci_size + +!===================================================================== +! LIST_KVII_SIZE counts the number of items in a list of (key,value) +! pairs +!===================================================================== + + function list_kvii_size(list) result(n) + + type(ListKeyValueII), pointer :: list + type(ListKeyValueII), pointer :: current + integer :: n + + if (associated(list)) then + n = 1 + current => list + do while (associated(current%next)) + current => current%next + n = n + 1 + enddo + else + n = 0 + endif + + end function list_kvii_size + +!===================================================================== +! LIST_INT_INSERT inserts a new element ! Arguments: ! elem Element in the linked list after ! which to insert the new element ! data The data for the new element !===================================================================== - subroutine list_insert( elem, data ) + subroutine list_int_insert(elem, data) - type(LinkedList), pointer :: elem - type(ListData), intent(in) :: data + type(ListInt), pointer :: elem + integer, intent(in) :: data - type(LinkedList), pointer :: next + type(ListInt), pointer :: next allocate(next) @@ -153,21 +336,90 @@ contains elem%next => next next%data = data - end subroutine list_insert + end subroutine list_int_insert !===================================================================== -! LIST_INSERT_HEAD inserts a new element before the first element +! LIST_REAL_INSERT inserts a new element +! Arguments: +! elem Element in the linked list after +! which to insert the new element +! data The data for the new element +!===================================================================== + + subroutine list_real_insert(elem, data) + + type(ListReal), pointer :: elem + real(8), intent(in) :: data + + type(ListReal), pointer :: next + + allocate(next) + + next%next => elem%next + elem%next => next + next%data = data + + end subroutine list_real_insert + +!===================================================================== +! LIST_KVCI_INSERT inserts a new element +! Arguments: +! elem Element in the linked list after +! which to insert the new element +! data The data for the new element +!===================================================================== + + subroutine list_kvci_insert(elem, data) + + type(ListKeyValueCI), pointer :: elem + type(KeyValueCI), intent(in) :: data + + type(ListKeyValueCI), pointer :: next + + allocate(next) + + next%next => elem%next + elem%next => next + next%data = data + + end subroutine list_kvci_insert + +!===================================================================== +! LIST_KVII_INSERT inserts a new element +! Arguments: +! elem Element in the linked list after +! which to insert the new element +! data The data for the new element +!===================================================================== + + subroutine list_kvii_insert(elem, data) + + type(ListKeyValueII), pointer :: elem + type(KeyValueII), intent(in) :: data + + type(ListKeyValueII), pointer :: next + + allocate(next) + + next%next => elem%next + elem%next => next + next%data = data + + end subroutine list_kvii_insert + +!===================================================================== +! LIST_INT_INSERT_HEAD inserts a new element before the first element ! Arguments: ! list Start of the list ! data The data for the new element !===================================================================== - subroutine list_insert_head( list, data ) + subroutine list_int_insert_head(list, data) - type(LinkedList), pointer :: list - type(ListData), intent(in) :: data + type(ListInt), pointer :: list + integer, intent(in) :: data - type(LinkedList), pointer :: elem + type(ListInt), pointer :: elem allocate(elem) elem%data = data @@ -175,33 +427,99 @@ contains elem%next => list list => elem - end subroutine list_insert_head + end subroutine list_int_insert_head !===================================================================== -! LIST_DELETE_ELEMENT deletes an element from the list +! LIST_REAL_INSERT_HEAD inserts a new element before the first element +! Arguments: +! list Start of the list +! data The data for the new element +!===================================================================== + + subroutine list_real_insert_head(list, data) + + type(ListReal), pointer :: list + real(8), intent(in) :: data + + type(ListReal), pointer :: elem + + allocate(elem) + elem%data = data + + elem%next => list + list => elem + + end subroutine list_real_insert_head + +!===================================================================== +! LIST_KVCI_INSERT_HEAD inserts a new element before the first element +! Arguments: +! list Start of the list +! data The data for the new element +!===================================================================== + + subroutine list_kvci_insert_head(list, data) + + type(ListKeyValueCI), pointer :: list + type(KeyValueCI), intent(in) :: data + + type(ListKeyValueCI), pointer :: elem + + allocate(elem) + elem%data = data + + elem%next => list + list => elem + + end subroutine list_kvci_insert_head + +!===================================================================== +! LIST_KVII_INSERT_HEAD inserts a new element before the first element +! Arguments: +! list Start of the list +! data The data for the new element +!===================================================================== + + subroutine list_kvii_insert_head(list, data) + + type(ListKeyValueII), pointer :: list + type(KeyValueII), intent(in) :: data + + type(ListKeyValueII), pointer :: elem + + allocate(elem) + elem%data = data + + elem%next => list + list => elem + + end subroutine list_kvii_insert_head + +!===================================================================== +! LIST_INT_DEL_ELEMENT deletes an element from the list ! Arguments: ! list Header of the list ! elem Element in the linked list to be removed !===================================================================== - subroutine list_delete_element( list, elem ) + subroutine list_int_del_element(list, elem) - type(LinkedList), pointer :: list - type(LinkedList), pointer :: elem + type(ListInt), pointer :: list + type(ListInt), pointer :: elem - type(LinkedList), pointer :: current - type(LinkedList), pointer :: prev + type(ListInt), pointer :: current + type(ListInt), pointer :: prev - if ( associated(list,elem) ) then + if (associated(list,elem)) then list => elem%next - deallocate( elem ) + deallocate(elem) else current => list prev => list - do while ( associated(current) ) - if ( associated(current,elem) ) then + do while (associated(current)) + if (associated(current,elem)) then prev%next => current%next - deallocate( current ) ! Is also "elem" + deallocate(current) ! Is also "elem" exit endif prev => current @@ -209,106 +527,190 @@ contains enddo endif - end subroutine list_delete_element + end subroutine list_int_del_element !===================================================================== -! LIST_GET_DATA gets the data stored with a list element +! LIST_REAL_DEL_ELEMENT deletes an element from the list ! Arguments: -! elem Element in the linked list +! list Header of the list +! elem Element in the linked list to be removed !===================================================================== - function list_get_data( elem ) result(data) + subroutine list_real_del_element(list, elem) - type(LinkedList), pointer :: elem - type(ListData) :: data + type(ListReal), pointer :: list + type(ListReal), pointer :: elem - data = elem%data + type(ListReal), pointer :: current + type(ListReal), pointer :: prev - end function list_get_data + if (associated(list,elem)) then + list => elem%next + deallocate(elem) + else + current => list + prev => list + do while (associated(current)) + if (associated(current,elem)) then + prev%next => current%next + deallocate(current) ! Is also "elem" + exit + endif + prev => current + current => current%next + enddo + endif + + end subroutine list_real_del_element !===================================================================== -! LIST_PUT_DATA stores new data with a list element +! LIST_KVCI_DEL_ELEMENT deletes an element from the list ! Arguments: -! elem Element in the linked list -! data The data to be stored +! list Header of the list +! elem Element in the linked list to be removed !===================================================================== - subroutine list_put_data( elem, data ) + subroutine list_kvci_del_element(list, elem) - type(LinkedList), pointer :: elem - type(ListData), intent(in) :: data + type(ListKeyValueCI), pointer :: list + type(ListKeyValueCI), pointer :: elem - elem%data = data + type(ListKeyValueCI), pointer :: current + type(ListKeyValueCI), pointer :: prev - end subroutine list_put_data + if (associated(list,elem)) then + list => elem%next + deallocate(elem) + else + current => list + prev => list + do while (associated(current)) + if (associated(current,elem)) then + prev%next => current%next + deallocate(current) ! Is also "elem" + exit + endif + prev => current + current => current%next + enddo + endif + + end subroutine list_kvci_del_element !===================================================================== -! DICT_CREATE creates and initializes a dictionary +! LIST_KVII_DEL_ELEMENT deletes an element from the list ! Arguments: -! dict Pointer to new dictionary -! key Key for the first element -! value Value for the first element -! Note: -! This version assumes a shallow copy is enough -! (that is, there are no pointers within the data -! to be stored) -! It also assumes the argument list does not already -! refer to a list. Use dict_destroy first to -! destroy up an old list. +! list Header of the list +! elem Element in the linked list to be removed !===================================================================== - subroutine dict_create( dict ) !, key, value ) + subroutine list_kvii_del_element(list, elem) - type(Dictionary), pointer :: dict -!!$ character(len=*), intent(in) :: key -!!$ integer, intent(in) :: value + type(ListKeyValueII), pointer :: list + type(ListKeyValueII), pointer :: elem -!!$ type(ListData) :: data - integer :: i -!!$ integer :: hash + type(ListKeyValueII), pointer :: current + type(ListKeyValueII), pointer :: prev - allocate( dict ) - allocate( dict%table(hash_size) ) + if (associated(list,elem)) then + list => elem%next + deallocate(elem) + else + current => list + prev => list + do while (associated(current)) + if (associated(current,elem)) then + prev%next => current%next + deallocate(current) ! Is also "elem" + exit + endif + prev => current + current => current%next + enddo + endif + + end subroutine list_kvii_del_element + +!===================================================================== +! DICT_CI_CREATE creates and initializes a dictionary +!===================================================================== + + subroutine dict_ci_create(dict) + + type(DictionaryCI), pointer :: dict + + integer :: i + + allocate(dict) + allocate(dict%table(hash_size)) do i = 1,hash_size dict%table(i)%list => null() enddo -!!$ data%key = key -!!$ data%value = value -!!$ -!!$ hash = dict_hashkey( trim(key ) ) -!!$ call list_create( dict%table(hash)%list, data ) - - end subroutine dict_create + end subroutine dict_ci_create !===================================================================== -! DICT_DESTROY destroys an entire dictionary -! Arguments: -! dict Pointer to the dictionary to be destroyed -! Note: -! This version assumes that there are no -! pointers within the data that need deallocation +! DICT_II_CREATE creates and initializes a dictionary !===================================================================== - subroutine dict_destroy( dict ) + subroutine dict_ii_create(dict) - type(Dictionary), pointer :: dict + type(DictionaryII), pointer :: dict - integer :: i + integer :: i + + allocate(dict) + allocate(dict%table(hash_size)) + + do i = 1,hash_size + dict%table(i)%list => null() + enddo + + end subroutine dict_ii_create + +!===================================================================== +! DICT_CI_DELETE destroys an entire dictionary +!===================================================================== + + subroutine dict_ci_delete(dict) + + type(DictionaryCI), pointer :: dict + + integer :: i do i = 1,size(dict%table) - if ( associated( dict%table(i)%list ) ) then - call list_destroy( dict%table(i)%list ) + if (associated(dict%table(i)%list)) then + call list_delete(dict%table(i)%list) endif enddo - deallocate( dict%table ) - deallocate( dict ) + deallocate(dict%table) + deallocate(dict) - end subroutine dict_destroy + end subroutine dict_ci_delete !===================================================================== -! DICT_ADD_KEY add a new keys +! DICT_II_DELETE destroys an entire dictionary +!===================================================================== + + subroutine dict_ii_delete(dict) + + type(DictionaryII), pointer :: dict + + integer :: i + + do i = 1,size(dict%table) + if (associated(dict%table(i)%list)) then + call list_delete(dict%table(i)%list) + endif + enddo + deallocate(dict%table) + deallocate(dict) + + end subroutine dict_ii_delete + +!===================================================================== +! DICT_CI_ADD_KEY add a new keys ! Arguments: ! dict Pointer to the dictionary ! key Key for the new element @@ -318,154 +720,313 @@ contains ! key's value is simply replaced !===================================================================== - subroutine dict_add_key( dict, key, value ) + subroutine dict_ci_add_key(dict, key, value) - type(Dictionary), pointer :: dict - character(len=*), intent(in) :: key - integer, intent(in) :: value + type(DictionaryCI), pointer :: dict + character(len=*), intent(in) :: key + integer, intent(in) :: value - type(ListData) :: data - type(LinkedList), pointer :: elem - integer :: hash + type(KeyValueCI) :: data + type(ListKeyValueCI), pointer :: elem + integer :: hash - elem => dict_get_elem( dict, key ) + elem => dict_get_elem(dict, key) - if ( associated(elem) ) then + if (associated(elem)) then elem%data%value = value else data%key = key data%value = value - hash = dict_hashkey( trim(key) ) - if ( associated( dict%table(hash)%list ) ) then - call list_insert( dict%table(hash)%list, data ) + hash = dict_hashkey(trim(key)) + if (associated(dict%table(hash)%list)) then + call list_insert(dict%table(hash)%list, data) else - call list_create( dict%table(hash)%list, data ) + call list_create(dict%table(hash)%list, data) endif endif - end subroutine dict_add_key + end subroutine dict_ci_add_key !===================================================================== -! DICT_DELETE_KEY deletes a key-value pair from the dictionary +! DICT_II_ADD_KEY add a new keys +! Arguments: +! dict Pointer to the dictionary +! key Key for the new element +! value Value for the new element +! Note: +! If the key already exists, the +! key's value is simply replaced +!===================================================================== + + subroutine dict_ii_add_key(dict, key, value) + + type(DictionaryII), pointer :: dict + integer, intent(in) :: key + integer, intent(in) :: value + + type(KeyValueII) :: data + type(ListKeyValueII), pointer :: elem + integer :: hash + + elem => dict_get_elem(dict, key) + + if (associated(elem)) then + elem%data%value = value + else + data%key = key + data%value = value + hash = dict_hashkey(key) + if (associated(dict%table(hash)%list)) then + call list_insert(dict%table(hash)%list, data) + else + call list_create(dict%table(hash)%list, data) + endif + endif + + end subroutine dict_ii_add_key + +!===================================================================== +! DICT_CI_DELETE_KEY deletes a key-value pair from the dictionary ! Arguments: ! dict Dictionary in question ! key Key to be removed !===================================================================== - subroutine dict_delete_key( dict, key ) + subroutine dict_ci_delete_key(dict, key) - type(Dictionary), pointer :: dict + type(DictionaryCI), pointer :: dict character(len=*), intent(in) :: key - type(LinkedList), pointer :: elem - integer :: hash + type(ListKeyValueCI), pointer :: elem + integer :: hash - elem => dict_get_elem( dict, key ) + elem => dict_get_elem(dict, key) - if ( associated(elem) ) then - hash = dict_hashkey( trim(key) ) - call list_delete_element( dict%table(hash)%list, elem ) + if (associated(elem)) then + hash = dict_hashkey(trim(key)) + call list_delete_element(dict%table(hash)%list, elem) endif - end subroutine dict_delete_key + end subroutine dict_ci_delete_key !===================================================================== -! DICT_GET_KEY gets the value belonging to a key +! DICT_II_DELETE_KEY deletes a key-value pair from the dictionary +! Arguments: +! dict Dictionary in question +! key Key to be removed +!===================================================================== + + subroutine dict_ii_delete_key(dict, key) + + type(DictionaryII), pointer :: dict + integer, intent(in) :: key + + type(ListKeyValueII), pointer :: elem + integer :: hash + + elem => dict_get_elem(dict, key) + + if (associated(elem)) then + hash = dict_hashkey(key) + call list_delete_element(dict%table(hash)%list, elem) + endif + + end subroutine dict_ii_delete_key + +!===================================================================== +! DICT_CI_GET_KEY gets the value belonging to a key ! Arguments: ! dict Pointer to the dictionary ! key Key for which the values are sought !===================================================================== - function dict_get_key( dict, key ) result(value) + function dict_ci_get_key(dict, key) result(value) - type(Dictionary), pointer :: dict - character(len=*), intent(in) :: key - integer :: value + type(DictionaryCI), pointer :: dict + character(len=*), intent(in) :: key + integer :: value - type(ListData) :: data - type(LinkedList), pointer :: elem + type(KeyValueCI) :: data + type(ListKeyValueCI), pointer :: elem - elem => dict_get_elem( dict, key ) + elem => dict_get_elem(dict, key) - if ( associated(elem) ) then + if (associated(elem)) then value = elem%data%value else value = DICT_NULL endif - end function dict_get_key + end function dict_ci_get_key !===================================================================== -! DICT_HAS_KEY checks if the dictionary has a particular key +! DICT_II_GET_KEY gets the value belonging to a key +! Arguments: +! dict Pointer to the dictionary +! key Key for which the values are sought +!===================================================================== + + function dict_ii_get_key(dict, key) result(value) + + type(DictionaryII), pointer :: dict + integer, intent(in) :: key + integer :: value + + type(KeyValueII) :: data + type(ListKeyValueII), pointer :: elem + + elem => dict_get_elem(dict, key) + + if (associated(elem)) then + value = elem%data%value + else + value = DICT_NULL + endif + + end function dict_ii_get_key + +!===================================================================== +! DICT_CI_HAS_KEY checks if the dictionary has a particular key ! Arguments: ! dict Pointer to the dictionary ! key Key to be sought !===================================================================== - function dict_has_key( dict, key ) result(has) + function dict_ci_has_key(dict, key) result(has) - type(Dictionary), pointer :: dict - character(len=*), intent(in) :: key - logical :: has + type(DictionaryCI), pointer :: dict + character(len=*), intent(in) :: key + logical :: has - type(LinkedList), pointer :: elem + type(ListKeyValueCI), pointer :: elem - elem => dict_get_elem( dict, key ) + elem => dict_get_elem(dict, key) has = associated(elem) - end function dict_has_key + end function dict_ci_has_key !===================================================================== -! DICT_GET_ELEM finds the element with a particular key +! DICT_II_HAS_KEY checks if the dictionary has a particular key ! Arguments: ! dict Pointer to the dictionary ! key Key to be sought !===================================================================== - function dict_get_elem( dict, key ) result(elem) + function dict_ii_has_key(dict, key) result(has) - type(Dictionary), pointer :: dict - character(len=*), intent(in) :: key + type(DictionaryII), pointer :: dict + integer, intent(in) :: key + logical :: has - type(LinkedList), pointer :: elem - integer :: hash + type(ListKeyValueII), pointer :: elem - hash = dict_hashkey( trim(key) ) + elem => dict_get_elem(dict, key) + + has = associated(elem) + + end function dict_ii_has_key + +!===================================================================== +! DICT_CI_GET_ELEM finds the element with a particular key +! Arguments: +! dict Pointer to the dictionary +! key Key to be sought +!===================================================================== + + function dict_ci_get_elem(dict, key) result(elem) + + type(DictionaryCI), pointer :: dict + character(len=*), intent(in) :: key + + type(ListKeyValueCI), pointer :: elem + integer :: hash + + hash = dict_hashkey(trim(key)) elem => dict%table(hash)%list - do while ( associated(elem) ) - if ( elem%data%key .eq. key ) then + do while (associated(elem)) + if (elem%data%key .eq. key) then exit else - elem => list_next( elem ) + elem => elem%next endif enddo - end function dict_get_elem + end function dict_ci_get_elem !===================================================================== -! DICT_HASHKEY determine the hash value from the string +! DICT_II_GET_ELEM finds the element with a particular key +! Arguments: +! dict Pointer to the dictionary +! key Key to be sought +!===================================================================== + + function dict_ii_get_elem(dict, key) result(elem) + + type(DictionaryII), pointer :: dict + integer, intent(in) :: key + + type(ListKeyValueII), pointer :: elem + integer :: hash + + hash = dict_hashkey(key) + elem => dict%table(hash)%list + do while (associated(elem)) + if (elem%data%key .eq. key) then + exit + else + elem => elem%next + endif + enddo + + end function dict_ii_get_elem + +!===================================================================== +! DICT_CI_HASHKEY determine the hash value from the string ! Arguments: ! key String to be examined !===================================================================== - integer function dict_hashkey( key ) + function dict_ci_hashkey(key) result(val) character(len=*), intent(in) :: key - integer :: hash - integer :: i + integer :: val + integer :: hash + integer :: i - dict_hashkey = 0 + val = 0 do i = 1,len(key) - dict_hashkey = multiplier * dict_hashkey + ichar(key(i:i)) + val = multiplier * val + ichar(key(i:i)) enddo - ! Added the absolute value on dict_hashkey-1 since the sum in the - ! do loop is susceptible to integer overflow - dict_hashkey = 1 + mod( abs(dict_hashkey-1), hash_size ) + ! Added the absolute val on val-1 since the sum in the do loop is + ! susceptible to integer overflow + val = 1 + mod(abs(val-1), hash_size) - end function dict_hashkey + end function dict_ci_hashkey + +!===================================================================== +! DICT_II_HASHKEY determine the hash value from the string +! Arguments: +! key String to be examined +!===================================================================== + + function dict_ii_hashkey(key) result(val) + + integer, intent(in) :: key + + integer :: val + integer :: hash + integer :: i + + val = 0 + + ! Added the absolute val on val-1 since the sum in the do loop is + ! susceptible to integer overflow + val = 1 + mod(abs(key-1), hash_size) + + end function dict_ii_hashkey end module data_structures diff --git a/src/energy_grid.f90 b/src/energy_grid.f90 index 6357e4bdf8..212827252a 100644 --- a/src/energy_grid.f90 +++ b/src/energy_grid.f90 @@ -2,6 +2,8 @@ module energy_grid use global use output, only: message + use data_structures, only: list_insert, list_size, list_delete, & + & list_size contains @@ -15,8 +17,8 @@ contains subroutine unionized_grid() - type(LinkedListGrid), pointer :: list => null() - type(LinkedListGrid), pointer :: current => null() + type(ListReal), pointer :: list => null() + type(ListReal), pointer :: current => null() type(Material), pointer :: mat => null() type(AceContinuous), pointer :: table => null() @@ -25,7 +27,7 @@ contains integer :: n character(100) :: msg - msg = "Creating unionized energy grid" + msg = "Creating unionized energy grid..." call message(msg, 5) ! loop over all materials @@ -42,6 +44,18 @@ contains end do end do + ! create allocated array from linked list + n_grid = list_size(list) + allocate(e_grid(n_grid)) + current => list + do i = 1,n_grid + e_grid(i) = current%data + current => current%next + end do + + ! delete linked list + call list_delete(list) + end subroutine unionized_grid !===================================================================== @@ -51,14 +65,14 @@ contains subroutine add_grid_points(list, energy, n_energy) - type(LinkedListGrid), pointer :: list + type(ListReal), pointer :: list real(8), intent(in) :: energy(n_energy) integer, intent(in) :: n_energy - type(LinkedListGrid), pointer :: current => null() - type(LinkedListGrid), pointer :: previous => null() - type(LinkedListGrid), pointer :: head => null() - type(LinkedListGrid), pointer :: tmp => null() + type(ListReal), pointer :: current => null() + type(ListReal), pointer :: previous => null() + type(ListReal), pointer :: head => null() + type(ListReal), pointer :: tmp => null() integer :: index real(8) :: E @@ -66,11 +80,11 @@ contains ! if the original list is empty, we need to allocate the first ! element and store first energy point - if (grid_count(list) == 0) then + if (list_size(list) == 0) then allocate(list) current => list do index = 1, n_energy - current%energy = energy(index) + current%data = energy(index) if (index == n_energy) then current%next => null() return @@ -93,7 +107,7 @@ contains do while (index <= n_energy) allocate(previous%next) current => previous%next - current%energy = energy(index) + current%data = energy(index) previous => current index = index + 1 end do @@ -101,10 +115,10 @@ contains exit end if - if (E < current%energy) then + if (E < current%data) then ! create new element and insert it in energy grid list allocate(tmp) - tmp%energy = E + tmp%data = E tmp%next => current if (associated(previous)) then previous%next => tmp @@ -119,7 +133,7 @@ contains index = index + 1 E = energy(index) - elseif (E == current%energy) then + elseif (E == current%data) then ! found the exact same energy, no need to store duplicates ! so just skip and move to next index index = index + 1 @@ -139,29 +153,37 @@ contains end subroutine add_grid_points !===================================================================== -! GRID_COUNT gives the number of energy points in list. This should -! eventually be replaced with a generic list function. +! ORIGINAL_INDICES !===================================================================== - function grid_count(list) result(count) + subroutine original_indices() - type(LinkedListGrid), pointer :: list - integer :: count + integer :: i, j + integer :: index + type(AceContinuous), pointer :: acecont - type(LinkedListGrid), pointer :: current + real(8) :: union_energy + real(8) :: energy + - ! determine size of list - if (associated(list)) then - count = 1 - current => list - do while (associated(current%next)) - current => current%next - count = count + 1 - enddo - else - count = 0 - end if + do i = 1, n_continuous + acecont => xs_continuous(i) + allocate(acecont%grid_index(n_grid)) - end function grid_count + index = 1 + energy = acecont%energy(index) + + do j = 1, n_grid + union_energy = e_grid(j) + if (union_energy >= energy) then + index = index + 1 + energy = acecont%energy(index) + end if + acecont%grid_index(j) = index-1 + end do + + end do + + end subroutine original_indices end module energy_grid diff --git a/src/fileio.f90 b/src/fileio.f90 index 8ae87126b2..59f520bc1a 100644 --- a/src/fileio.f90 +++ b/src/fileio.f90 @@ -2,9 +2,11 @@ module fileio use global use types, only: Cell, Surface, ExtSource - use string, only: split_string_wl, lower_case, str_to_real, split_string + use string, only: split_string_wl, lower_case, str_to_real, split_string, & + & concatenate use output, only: message, warning, error - use data_structures, only: Dictionary, dict_create, dict_add_key, dict_get_key + use data_structures, only: dict_create, dict_add_key, dict_get_key, & + & DICT_NULL implicit none @@ -45,7 +47,10 @@ contains if (.not. file_exists) then msg = "Input file '" // trim(path_input) // "' does not exist!" call error(msg) - elseif (readable(1:3) /= 'YES') then + elseif (readable(1:3) == 'NO') then + ! Need to explicitly check for a NO status -- Intel compiler + ! looks at file attributes only if the file is open on a + ! unit. Therefore, it will always return UNKNOWN msg = "Input file '" // trim(path_input) // "' is not readable! & &Change file permissions with chmod command." call error(msg) @@ -66,21 +71,26 @@ contains character(250) :: line, msg character(32) :: words(max_words) integer :: in = 7 - integer :: readError + integer :: ioError integer :: n + msg = "First pass through input file..." + call message(msg, 5) + n_cells = 0 n_surfaces = 0 n_materials = 0 - open(file=filename, unit=in, status='old', action='read') + open(FILE=filename, UNIT=in, STATUS='old', & + & ACTION='read', IOSTAT=ioError) + if (ioError /= 0) then + msg = "Error while opening file: " // filename + call error(msg) + end if do - read(unit=in, fmt='(A250)', iostat=readError) line - if (readError /= 0) exit - ! TODO: don't really need to split entire string, this could be - ! replaced by just finding first word - call split_string(line, words, n) + call get_next_line(in, words, n, ioError) + if (ioError /= 0) exit if (n == 0) cycle select case (trim(words(1))) @@ -97,19 +107,19 @@ contains ! defined for the problem if (n_cells == 0) then msg = "No cells specified!" - close(unit=in) + close(UNIT=in) call error(msg) elseif (n_surfaces == 0) then msg = "No surfaces specified!" - close(unit=in) + close(UNIT=in) call error(msg) elseif (n_materials == 0) then msg = "No materials specified!" - close(unit=in) + close(UNIT=in) call error(msg) end if - close(unit=in) + close(UNIT=in) ! Allocate arrays for cells, surfaces, etc. allocate(cells(n_cells)) @@ -119,6 +129,7 @@ contains ! Set up dictionaries call dict_create(cell_dict) call dict_create(surface_dict) + call dict_create(material_dict) end subroutine read_count @@ -132,7 +143,7 @@ contains character(*), intent(in) :: filename - character(250) :: line + character(250) :: line, msg character(32) :: words(max_words) integer :: in = 7 integer :: ioError @@ -142,17 +153,24 @@ contains integer :: index_material integer :: index_source + msg = "Second pass through input file..." + call message(msg, 5) + index_cell = 0 index_surface = 0 index_material = 0 index_source = 0 - open(file=filename, unit=in, status='old', action='read') + open(FILE=filename, UNIT=in, STATUS='old', & + & ACTION='read', IOSTAT=ioError) + if (ioError /= 0) then + msg = "Error while opening file: " // filename + call error(msg) + end if do - read(unit=in, fmt='(A250)', iostat=ioError) line + call get_next_line(in, words, n, ioError) if ( ioError /= 0 ) exit - call split_string_wl(line, words, n) ! skip comments if (n==0 .or. words(1)(1:1) == '#') cycle @@ -187,12 +205,57 @@ contains end do - close(unit=in) + close(UNIT=in) - ! call adjust_cell_index + call adjust_indices() end subroutine read_input +!===================================================================== +! ADJUST_INDICES changes the boundary_list values for each cell and +! the material index assigned to each to the indices in the surfaces +! and material array rather than the unique IDs assigned to each +! surface and material +!===================================================================== + + subroutine adjust_indices() + + type(Cell), pointer :: c + + integer :: i, j + integer :: index + integer :: surf_num + character(250) :: msg + + do i = 1, n_cells + ! adjust boundary list + c => cells(i) + do j = 1, c%n_items + surf_num = c%boundary_list(j) + if (surf_num < OP_DIFFERENCE) then + index = dict_get_key(surface_dict, abs(surf_num)) + if (index == DICT_NULL) then + surf_num = abs(surf_num) + msg = "Could not find surface " // trim(int_to_str(surf_num)) // & + & " specified on cell " // trim(int_to_str(c%uid)) + call error(msg) + end if + c%boundary_list(j) = sign(index, surf_num) + end if + end do + + ! adjust material indices + index = dict_get_key(material_dict, c%material) + if (index == DICT_NULL) then + msg = "Could not find material " // trim(int_to_str(c%material)) // & + & " specified on cell " // trim(int_to_str(c%uid)) + call error(msg) + end if + c%material = index + end do + + end subroutine adjust_indices + !===================================================================== ! READ_CELL parses the data on a cell card. !===================================================================== @@ -203,8 +266,9 @@ contains character(*), intent(in) :: words(max_words) integer, intent(in) :: n_words - integer :: readError + integer :: ioError integer :: i + integer :: n_items character(250) :: msg character(32) :: word type(Cell), pointer :: this_cell => null() @@ -212,19 +276,23 @@ contains this_cell => cells(index) ! Read cell identifier - read(words(2), fmt='(I8)', iostat=readError) this_cell%uid - if (readError > 0) then + read(words(2), FMT='(I8)', IOSTAT=ioError) this_cell%uid + if (ioError > 0) then msg = "Invalid cell name: " // words(2) call error(msg) end if - call dict_add_key(cell_dict, words(2), index) + call dict_add_key(cell_dict, this_cell%uid, index) ! Read cell material - read(words(3), fmt='(I8)', iostat=readError) this_cell%material + read(words(3), FMT='(I8)', IOSTAT=ioError) this_cell%material + + ! Assign number of items + n_items = n_words - 3 + this_cell%n_items = n_items ! Read list of surfaces - allocate(this_cell%boundary_list(n_words-3)) - do i = 1, n_words-3 + allocate(this_cell%boundary_list(n_items)) + do i = 1, n_items word = words(i+3) if (word(1:1) == '(') then this_cell%boundary_list(i) = OP_LEFT_PAREN @@ -235,7 +303,7 @@ contains elseif (word(1:1) == '#') then this_cell%boundary_list(i) = OP_DIFFERENCE else - read(word, fmt='(I8)', iostat=readError) this_cell%boundary_list(i) + read(word, FMT='(I8)', IOSTAT=ioError) this_cell%boundary_list(i) end if end do @@ -251,7 +319,7 @@ contains character(*), intent(in) :: words(max_words) integer, intent(in) :: n_words - integer :: readError + integer :: ioError integer :: i integer :: coeffs_reqd character(250) :: msg @@ -261,11 +329,12 @@ contains this_surface => surfaces(index) ! Read surface identifier - read(words(2), fmt='(I8)', iostat=readError) this_surface%uid - if (readError > 0) then + read(words(2), FMT='(I8)', IOSTAT=ioError) this_surface%uid + if (ioError > 0) then msg = "Invalid surface name: " // words(2) call error(msg) end if + call dict_add_key(surface_dict, this_surface%uid, index) ! Read surface type word = words(3) @@ -340,7 +409,7 @@ contains character(250) :: msg character(32) :: word - integer :: readError + integer :: ioError integer :: values_reqd integer :: i @@ -371,8 +440,8 @@ contains end subroutine read_source !===================================================================== -! READ_MATERIAL parses a material card, creating the appropriate -! s and s. +! READ_MATERIAL parses a material card. Note that atom percents and +! densities are normalized in a separate routine !===================================================================== subroutine read_material(index, words, n_words) @@ -383,18 +452,30 @@ contains character(100) :: msg integer :: i + integer :: ioError integer :: n_isotopes type(Material), pointer :: mat - if (mod(n_words,2) /= 0) then + if (mod(n_words,2) == 0 .or. n_words < 5) then msg = "Invalid number of arguments for material: " // words(2) call error(msg) end if - n_isotopes = (n_words-2)/2 + n_isotopes = (n_words-3)/2 mat => materials(index) mat%n_isotopes = n_isotopes + ! Read surface identifier + read(words(2), FMT='(I8)', IOSTAT=ioError) mat%uid + if (ioError > 0) then + msg = "Invalid surface name: " // words(2) + call error(msg) + end if + call dict_add_key(material_dict, mat%uid, index) + + ! Read atom density + mat%atom_density = str_to_real(words(3)) + ! allocate isotope and density list allocate(mat%names(n_isotopes)) allocate(mat%isotopes(n_isotopes)) @@ -403,8 +484,8 @@ contains ! read names and percentage do i = 1, n_isotopes - mat%names(i) = words(2*i+1) - mat%atom_percent(i) = str_to_real(words(2*i+2)) + mat%names(i) = words(2*i+2) + mat%atom_percent(i) = str_to_real(words(2*i+3)) end do end subroutine read_material @@ -420,9 +501,13 @@ contains type(Material), pointer :: mat integer :: index integer :: i, j - character(10) :: key real(8) :: sum_percent - real(8) :: awr, w + real(8) :: awr ! atomic weight ratio + real(8) :: w ! weight percent + real(8) :: x ! atom percent + logical :: percent_in_atom + logical :: density_in_atom + character(10) :: key character(100) :: msg ! first find the index in the xsdata array for each isotope in @@ -439,39 +524,48 @@ contains call error(msg) end if - ! If data given in a/o, normalize it by sum - if (mat%atom_percent(1) > 0) then - sum_percent = sum(mat%atom_percent) - mat%atom_percent = mat%atom_percent / sum_percent - do j = 1, mat%n_isotopes - key = mat%names(j) - index = dict_get_key(xsdata_dict, key) - mat%isotopes(j) = index - end do - cycle - end if + percent_in_atom = (mat%atom_percent(1) > 0) + density_in_atom = (mat%atom_density > 0) - ! Otherwise, need to find atomic weight ratios for each isotope - ! in order to normalize to a/o. To convert weight percent to - ! atom percent, the formula is: x = (w/awr) / sum(w/awr) sum_percent = 0 do j = 1, mat%n_isotopes + ! Set indices for isotopes key = mat%names(j) index = dict_get_key(xsdata_dict, key) mat%isotopes(j) = index - awr = xsdatas(index)%awr - w = mat%atom_percent(j) - sum_percent = sum_percent + w/awr - end do - - ! change each isotope from w/o to a/o - do j = 1, mat%n_isotopes - index = mat%isotopes(j) + ! determine atomic weight ratio awr = xsdatas(index)%awr - w = mat%atom_percent(j) - mat%atom_percent(j) = (w/awr)/sum_percent + + ! if given weight percent, convert all values so that they + ! are divided by awr. thus, when a sum is done over the + ! values, it's actually sum(w/awr) + if (.not. percent_in_atom) then + mat%atom_percent(j) = -mat%atom_percent(j) / awr + end if end do + + ! determine normalized atom percents. if given atom percents, + ! this is straightforward. if given weight percents, the value + ! is w/awr and is divided by sum(w/awr) + sum_percent = sum(mat%atom_percent) + mat%atom_percent = mat%atom_percent / sum_percent + + ! Change density in g/cm^3 to atom/b-cm. Since all values are + ! now in atom percent, the sum needs to be re-evaluated as + ! 1/sum(x*awr) + if (.not. density_in_atom) then + sum_percent = 0 + do j = 1, mat%n_isotopes + index = mat%isotopes(j) + awr = xsdatas(index)%awr + x = mat%atom_percent(j) + sum_percent = sum_percent + x*awr + end do + sum_percent = 1.0_8 / sum_percent + mat%atom_density = -mat%atom_density * N_AVOGADRO & + & / MASS_NEUTRON * sum_percent + end if end do end subroutine normalize_ao @@ -501,28 +595,28 @@ contains integer, intent(in) :: n_words character(250) :: msg - integer :: readError + integer :: ioError ! Set problem type to criticality problem_type = PROB_CRITICALITY ! Read number of cycles - read(words(2), fmt='(I8)', iostat=readError) n_cycles - if (readError > 0) then + read(words(2), FMT='(I8)', IOSTAT=ioError) n_cycles + if (ioError > 0) then msg = "Invalid number of cycles: " // words(2) call error(msg) end if ! Read number of inactive cycles - read(words(3), fmt='(I8)', iostat=readError) n_inactive - if (readError > 0) then + read(words(3), FMT='(I8)', IOSTAT=ioError) n_inactive + if (ioError > 0) then msg = "Invalid number of inactive cycles: " // words(2) call error(msg) end if ! Read number of particles - read(words(4), fmt='(I8)', iostat=readError) n_particles - if ( readError > 0 ) then + read(words(4), FMT='(I8)', IOSTAT=ioError) n_particles + if (ioError > 0) then msg = "Invalid number of particles: " // words(2) call error(msg) end if @@ -543,6 +637,43 @@ contains end subroutine read_line +!===================================================================== +! GET_NEXT_LINE reads the next line to the file connected on the +! specified unit including any continuation lines. If a line ends in +! an ampersand, the next line is read and its words are appended to +! the final array +!===================================================================== + + subroutine get_next_line(unit, words, n, ioError) + + integer, intent(in) :: unit + character(*), intent(out) :: words(max_words) + integer, intent(out) :: n + integer, intent(out) :: ioError + + character(250) :: line + character(32) :: local_words(max_words) + integer :: index + + index = 0 + do + read(UNIT=unit, FMT='(A250)', IOSTAT=ioError) line + if (ioError /= 0) return + call split_string_wl(line, local_words, n) + if (n == 0) exit + if (local_words(n) == '&') then + words(index+1:index+n-1) = local_words(1:n-1) + index = index + n-1 + else + words(index+1:index+n) = local_words(1:n) + index = index + n + exit + end if + end do + n = index + + end subroutine get_next_line + !===================================================================== ! SKIP_LINES skips 'n_lines' lines from a file open on a unit !===================================================================== @@ -557,7 +688,7 @@ contains character(max_line) :: tmp do i = 1, n_lines - read(UNIT=unit, fmt='(A250)', IOSTAT=ioError) tmp + read(UNIT=unit, FMT='(A250)', IOSTAT=ioError) tmp end do end subroutine skip_lines diff --git a/src/geometry.f90 b/src/geometry.f90 index 9f941cffc9..4aea337be4 100644 --- a/src/geometry.f90 +++ b/src/geometry.f90 @@ -3,6 +3,7 @@ module geometry use global use types, only: Cell, Surface use output, only: error, message + use data_structures, only: dict_get_key implicit none @@ -12,11 +13,11 @@ contains ! CELL_CONTAINS determines whether a given point is inside a cell !===================================================================== - function cell_contains( c, xyz ) + function cell_contains(c, neut) result(in_cell) - type(Cell), intent(in) :: c - real(8), intent(in) :: xyz(3) - logical :: cell_contains + type(Cell), pointer :: c + type(Neutron), pointer :: neut + logical :: in_cell integer, allocatable :: expression(:) integer :: specified_sense @@ -24,44 +25,42 @@ contains integer :: n_boundaries integer :: i, j integer :: surf_num + integer :: current_surface type(Surface), pointer :: surf => null() - logical :: surface_found character(250) :: msg + current_surface = neut%surface + n_boundaries = size(c%boundary_list) - allocate( expression(n_boundaries) ) + allocate(expression(n_boundaries)) expression = c%boundary_list do i = 1, n_boundaries - surface_found = .false. ! Don't change logical operator - if ( expression(i) >= OP_DIFFERENCE ) then + if (expression(i) >= OP_DIFFERENCE) then cycle end if ! Lookup surface - surf_num = abs(expression(i)) - ! TODO: replace this loop with a hash since this lookup is O(N) - do j = 1, n_surfaces - surf => surfaces(j) - if ( surf%uid == surf_num ) then - surface_found = .true. - exit - end if - end do + surf_num = expression(i) + surf => surfaces(abs(surf_num)) - ! Report error if can't find specified surface - if ( .not. surface_found ) then - deallocate( expression ) - msg = "Count not find surface " // trim(int_to_str(surf_num)) // & - & " specified on cell " // int_to_str(c%uid) - call error( msg ) + ! Check if the particle is currently on the specified surface + if (surf_num == current_surface) then + ! neutron is on specified surface heading into cell + expression(i) = 1 + cycle + elseif (surf_num == -current_surface) then + ! neutron is on specified surface, but heading other + ! direction + expression(i) = 0 + cycle end if ! Compare sense of point to specified sense specified_sense = sign(1,expression(i)) - call sense( surf, xyz, actual_sense ) - if ( actual_sense == specified_sense ) then + call sense(surf, neut%xyz, actual_sense) + if (actual_sense == specified_sense) then expression(i) = 1 else expression(i) = 0 @@ -71,14 +70,14 @@ contains ! TODO: Need to replace this with a 'lgeval' like subroutine that ! can actually test expressions with unions and parentheses - if ( all( expression == 1 ) ) then - cell_contains = .true. + if (all(expression == 1)) then + in_cell = .true. else - cell_contains = .false. + in_cell = .false. end if ! Free up memory from expression - deallocate( expression ) + deallocate(expression) end function cell_contains @@ -86,28 +85,34 @@ contains ! FIND_CELL determines what cell a source neutron is in !===================================================================== - subroutine find_cell( neut ) + subroutine find_cell(neut) type(Neutron), pointer, intent(inout) :: neut + type(Cell), pointer :: this_cell logical :: found_cell character(250) :: msg integer :: i ! determine what region in do i = 1, n_cells - if ( cell_contains(cells(i), neut%xyz) ) then + this_cell => cells(i) + if (cell_contains(this_cell, neut)) then neut%cell = i found_cell = .true. + + ! set current pointers + cCell => this_cell + cMaterial => materials(cCell%material) exit end if end do ! if neutron couldn't be located, print error - if ( .not. found_cell ) then + if (.not. found_cell) then write(msg, 100) "Could not locate cell for neutron at: ", neut%xyz 100 format (A,3ES11.3) - call error( msg ) + call error(msg) end if end subroutine find_cell @@ -116,45 +121,50 @@ contains ! CROSS_BOUNDARY moves a neutron into a new cell !===================================================================== - subroutine cross_boundary( neut ) + subroutine cross_boundary(neut) type(Neutron), pointer, intent(in) :: neut type(Surface), pointer :: surf - type(Cell), pointer :: c + type(Cell), pointer :: c integer :: i integer :: index_cell - character(100) :: msg + character(250) :: msg surf => surfaces(abs(neut%surface)) + if (verbosity >= 10) then + msg = " Crossing surface " // trim(int_to_str(neut%surface)) + call message(msg, 10) + end if ! check for leakage - if ( surf%bc == BC_VACUUM ) then + if (surf%bc == BC_VACUUM) then neut%alive = .false. - msg = "Particle " // trim(int_to_str(neut%uid)) // " leaked out of surface " & - & // trim(int_to_str(surf%uid)) - call message( msg, 10 ) + if (verbosity >= 10) then + msg = " Leaked out of surface " // trim(int_to_str(surf%uid)) + call message(msg, 10) + end if return end if - if ( neut%surface < 0 ) then + if (neut%surface < 0 .and. allocated(surf%neighbor_pos)) then ! If coming from negative side of surface, search all the ! neighboring cells on the positive side - do i = 1, size(surf%neighbor_pos(:)) + do i = 1, size(surf%neighbor_pos) index_cell = surf%neighbor_pos(i) - c = cells(index_cell) - if ( cell_contains(c, neut%xyz) ) then + c => cells(index_cell) + if (cell_contains(c, neut)) then neut%cell = index_cell return end if end do - elseif ( neut%surface > 0 ) then + elseif (neut%surface > 0 .and. allocated(surf%neighbor_neg)) then ! If coming from positive side of surface, search all the ! neighboring cells on the negative side - do i = 1, size(surf%neighbor_neg(:)) + do i = 1, size(surf%neighbor_neg) index_cell = surf%neighbor_neg(i) - c = cells(index_cell) - if ( cell_contains(c, neut%xyz) ) then + c => cells(index_cell) + if (cell_contains(c, neut)) then neut%cell = index_cell return end if @@ -164,17 +174,17 @@ contains ! Couldn't find particle in neighboring cells, search through all ! cells do i = 1, size(cells) - c = cells(i) - if ( cell_contains(c, neut%xyz) ) then + c => cells(i) + if (cell_contains(c, neut)) then neut%cell = i return end if end do ! Couldn't find next cell anywhere! - msg = "After neutron crossed surface " // int_to_str(neut%surface) // & + msg = "After neutron crossed surface " // trim(int_to_str(neut%surface)) // & & ", it could not be located in any cell and it did not leak." - call error( msg ) + call error(msg) end subroutine cross_boundary @@ -184,7 +194,7 @@ contains ! direction. !===================================================================== - subroutine dist_to_boundary( neut, dist, surf, other_cell ) + subroutine dist_to_boundary(neut, dist, surf, other_cell) type(Neutron), intent(in) :: neut real(8), intent(out) :: dist @@ -197,6 +207,7 @@ contains integer :: n_boundaries integer, allocatable :: expression(:) integer :: index_surf + integer :: current_surf real(8) :: x,y,z,u,v,w real(8) :: d real(8) :: x0,y0,z0,r @@ -205,12 +216,14 @@ contains real(8) :: quad character(250) :: msg - if ( present(other_cell) ) then + if (present(other_cell)) then cell_p => cells(other_cell) else cell_p => cells(neut%cell) end if + current_surf = neut%surface + x = neut%xyz(1) y = neut%xyz(2) z = neut%xyz(3) @@ -220,58 +233,63 @@ contains dist = INFINITY n_boundaries = size(cell_p%boundary_list) - allocate( expression(n_boundaries) ) + allocate(expression(n_boundaries)) expression = cell_p%boundary_list do i = 1, n_boundaries - index_surf = abs(expression(i)) - if ( index_surf >= OP_DIFFERENCE ) cycle + ! check for coincident surfaec + index_surf = expression(i) + if (index_surf == current_surf) cycle + + ! check for operators + index_surf = abs(index_surf) + if (index_surf >= OP_DIFFERENCE) cycle surf_p => surfaces(index_surf) - select case ( surf_p%type ) - case ( SURF_PX ) - if ( u == 0.0 ) then + select case (surf_p%type) + case (SURF_PX) + if (u == 0.0) then d = INFINITY else x0 = surf_p%coeffs(1) d = (x0 - x)/u - if ( d < 0 ) d = INFINITY + if (d < 0) d = INFINITY end if - case ( SURF_PY ) - if ( v == 0.0 ) then + case (SURF_PY) + if (v == 0.0) then d = INFINITY else y0 = surf_p%coeffs(1) d = (y0 - y)/v - if ( d < 0 ) d = INFINITY + if (d < 0) d = INFINITY end if - case ( SURF_PZ ) - if ( w == 0.0 ) then + case (SURF_PZ) + if (w == 0.0) then d = INFINITY else z0 = surf_p%coeffs(1) d = (z0 - z)/w - if ( d < 0.0 ) d = INFINITY + if (d < 0.0) d = INFINITY end if - case ( SURF_PLANE ) + case (SURF_PLANE) A = surf_p%coeffs(1) B = surf_p%coeffs(2) C = surf_p%coeffs(3) D = surf_p%coeffs(4) tmp = A*u + B*v + C*w - if ( tmp == 0.0 ) then + if (tmp == 0.0) then d = INFINITY else d = -(A*x + B*y + C*w - D)/tmp - if ( d < 0.0 ) d = INFINITY + if (d < 0.0) d = INFINITY end if - case ( SURF_CYL_X ) + case (SURF_CYL_X) a = 1.0 - u**2 ! v^2 + w^2 - if ( a == 0.0 ) then + if (a == 0.0) then d = INFINITY else y0 = surf_p%coeffs(1) @@ -284,7 +302,7 @@ contains c = y**2 + z**2 - r**2 quad = k**2 - a*c - if ( c < 0 ) then + if (c < 0) then ! particle is inside the cylinder, thus one distance ! must be negative and one must be positive. The ! positive distance will be the one with negative sign @@ -299,14 +317,14 @@ contains ! positive sign on sqrt(quad) d = -(k + sqrt(quad))/a - if ( d < 0 ) d = INFINITY + if (d < 0) d = INFINITY end if end if - case ( SURF_CYL_Y ) + case (SURF_CYL_Y) a = 1.0 - v**2 ! u^2 + w^2 - if ( a == 0.0 ) then + if (a == 0.0) then d = INFINITY else x0 = surf_p%coeffs(1) @@ -319,7 +337,7 @@ contains c = x**2 + z**2 - r**2 quad = k**2 - a*c - if ( c < 0 ) then + if (c < 0) then ! particle is inside the cylinder, thus one distance ! must be negative and one must be positive. The ! positive distance will be the one with negative sign @@ -334,14 +352,14 @@ contains ! positive sign on sqrt(quad) d = -(k + sqrt(quad))/a - if ( d < 0 ) d = INFINITY + if (d < 0) d = INFINITY end if end if - case ( SURF_CYL_Z ) + case (SURF_CYL_Z) a = 1.0 - w**2 ! u^2 + v^2 - if ( a == 0.0 ) then + if (a == 0.0) then d = INFINITY else x0 = surf_p%coeffs(1) @@ -354,12 +372,12 @@ contains c = x**2 + y**2 - r**2 quad = k**2 - a*c - if ( quad < 0 ) then + if (quad < 0) then ! no intersection with cylinder d = INFINITY - elseif ( c < 0 ) then + elseif (c < 0) then ! particle is inside the cylinder, thus one distance ! must be negative and one must be positive. The ! positive distance will be the one with negative sign @@ -374,12 +392,12 @@ contains ! positive sign on sqrt(quad) d = -(k + sqrt(quad))/a - if ( d < 0 ) d = INFINITY + if (d < 0) d = INFINITY end if end if - case ( SURF_SPHERE ) + case (SURF_SPHERE) x0 = surf_p%coeffs(1) y0 = surf_p%coeffs(2) z0 = surf_p%coeffs(3) @@ -392,12 +410,12 @@ contains c = x**2 + y**2 + z**2 - r**2 quad = k**2 - c - if ( quad < 0 ) then + if (quad < 0) then ! no intersection with sphere d = INFINITY - elseif ( c < 0 ) then + elseif (c < 0) then ! particle is inside the sphere, thus one distance ! must be negative and one must be positive. The ! positive distance will be the one with negative sign @@ -412,26 +430,26 @@ contains ! positive sign on sqrt(quad) d = -(k + sqrt(quad)) - if ( d < 0 ) d = INFINITY + if (d < 0) d = INFINITY end if - case ( SURF_GQ ) + case (SURF_GQ) msg = "Surface distance not yet implement for general quadratic." - call error( msg ) + call error(msg) end select ! Check is calculated distance is new minimum - if ( d < dist ) then + if (d < dist) then dist = d - surf = expression(i) + surf = -expression(i) end if end do ! deallocate expression - deallocate( expression ) + deallocate(expression) end subroutine dist_to_boundary @@ -441,7 +459,7 @@ contains ! a particular point is in. !===================================================================== - subroutine sense( surf, xyz, s ) + subroutine sense(surf, xyz, s) type(Surface), intent(in) :: surf real(8), intent(in) :: xyz(3) @@ -457,110 +475,110 @@ contains y = xyz(2) z = xyz(3) - select case ( surf%type ) - case ( SURF_PX ) + select case (surf%type) + case (SURF_PX) x0 = surf%coeffs(1) func = x - x0 - case ( SURF_PY ) + case (SURF_PY) y0 = surf%coeffs(1) func = y - y0 - case ( SURF_PZ ) + case (SURF_PZ) z0 = surf%coeffs(1) func = z - z0 - case ( SURF_PLANE ) + case (SURF_PLANE) A = surf%coeffs(1) B = surf%coeffs(2) C = surf%coeffs(3) D = surf%coeffs(4) func = A*x + B*y + C*z - D - case ( SURF_CYL_X ) + case (SURF_CYL_X) y0 = surf%coeffs(1) z0 = surf%coeffs(2) r = surf%coeffs(3) func = (y-y0)**2 + (z-z0)**2 - r**2 - case ( SURF_CYL_Y ) + case (SURF_CYL_Y) x0 = surf%coeffs(1) z0 = surf%coeffs(2) r = surf%coeffs(3) func = (x-x0)**2 + (z-z0)**2 - r**2 - case ( SURF_CYL_Z ) + case (SURF_CYL_Z) x0 = surf%coeffs(1) y0 = surf%coeffs(2) r = surf%coeffs(3) func = (x-x0)**2 + (y-y0)**2 - r**2 - case ( SURF_SPHERE ) + case (SURF_SPHERE) x0 = surf%coeffs(1) y0 = surf%coeffs(2) z0 = surf%coeffs(3) r = surf%coeffs(4) func = (x-x0)**2 + (y-y0)**2 + (z-z0)**2 - r**2 - case ( SURF_BOX_X ) + case (SURF_BOX_X) y0 = surf%coeffs(1) z0 = surf%coeffs(2) y1 = surf%coeffs(3) z1 = surf%coeffs(4) - if ( y >= y0 .and. y < y1 .and. z >= z0 .and. z < z1 ) then + if (y >= y0 .and. y < y1 .and. z >= z0 .and. z < z1) then s = SENSE_NEGATIVE else s = SENSE_POSITIVE end if return - case ( SURF_BOX_Y ) + case (SURF_BOX_Y) x0 = surf%coeffs(1) z0 = surf%coeffs(2) x1 = surf%coeffs(3) z1 = surf%coeffs(4) - if ( x >= x0 .and. x < x1 .and. z >= z0 .and. z < z1 ) then + if (x >= x0 .and. x < x1 .and. z >= z0 .and. z < z1) then s = SENSE_NEGATIVE else s = SENSE_POSITIVE end if return - case ( SURF_BOX_Z ) + case (SURF_BOX_Z) x0 = surf%coeffs(1) y0 = surf%coeffs(2) x1 = surf%coeffs(3) y1 = surf%coeffs(4) - if ( x >= x0 .and. x < x1 .and. y >= y0 .and. y < y1 ) then + if (x >= x0 .and. x < x1 .and. y >= y0 .and. y < y1) then s = SENSE_NEGATIVE else s = SENSE_POSITIVE end if return - case ( SURF_BOX ) + case (SURF_BOX) x0 = surf%coeffs(1) y0 = surf%coeffs(2) z0 = surf%coeffs(3) x1 = surf%coeffs(4) y1 = surf%coeffs(5) z1 = surf%coeffs(6) - if ( x >= x0 .and. x < x1 .and. y >= y0 .and. y < y1 .and. & - & z >= z0 .and. z < z1 ) then + if (x >= x0 .and. x < x1 .and. y >= y0 .and. y < y1 .and. & + & z >= z0 .and. z < z1) then s = SENSE_NEGATIVE else s = SENSE_POSITIVE end if return - case ( SURF_GQ ) + case (SURF_GQ) func = A*x**2 + B*y**2 + C*z**2 + D*x*y + E*y*z + F*x*z + G*x & & + H*y + I*z + J end select ! Check which side of surface the point is on - if ( func > 0 ) then + if (func > 0) then s = SENSE_POSITIVE else s = SENSE_NEGATIVE @@ -568,4 +586,83 @@ contains end subroutine sense +!===================================================================== +! NEIGHBOR_LISTS builds a list of neighboring cells to each surface to +! speed up searches when a cell boundary is crossed. +!===================================================================== + + subroutine neighbor_lists() + + type(Cell), pointer :: c + type(Surface), pointer :: surf + integer :: i, j + integer :: index + integer :: surf_num + logical :: positive + character(250) :: msg + + integer, allocatable :: count_positive(:) + integer, allocatable :: count_negative(:) + + msg = "Building neighboring cells lists for each surface..." + call message(msg, 4) + + allocate(count_positive(n_surfaces)) + allocate(count_negative(n_surfaces)) + count_positive = 0 + count_negative = 0 + + do i = 1, n_cells + c => cells(i) + + ! loop over each surface specification + do j = 1, c%n_items + index = c%boundary_list(j) + positive = (index > 0) + index = abs(index) + if (positive) then + count_positive(index) = count_positive(index) + 1 + else + count_negative(index) = count_negative(index) + 1 + end if + end do + end do + + ! allocate neighbor lists for each surface + do i = 1, n_surfaces + surf => surfaces(i) + if (count_positive(i) > 0) then + allocate(surf%neighbor_pos(count_positive(i))) + end if + if (count_negative(i) > 0) then + allocate(surf%neighbor_neg(count_negative(i))) + end if + end do + + count_positive = 0 + count_negative = 0 + + ! loop over all cells + do i = 1, n_cells + c => cells(i) + + ! loop over each surface specification + do j = 1, c%n_items + index = c%boundary_list(j) + positive = (index > 0) + index = abs(index) + + surf => surfaces(index) + if (positive) then + count_positive(index) = count_positive(index) + 1 + surf%neighbor_pos(count_positive(index)) = i + else + count_negative(index) = count_negative(index) + 1 + surf%neighbor_neg(count_negative(index)) = i + end if + end do + end do + + end subroutine neighbor_lists + end module geometry diff --git a/src/global.f90 b/src/global.f90 index 5b6c3ae431..ccc63fccc5 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -14,18 +14,28 @@ module global integer :: n_surfaces ! # of surfaces integer :: n_materials ! # of materials - type(Dictionary), pointer :: cell_dict - type(Dictionary), pointer :: surface_dict - type(Dictionary), pointer :: xsdata_dict - type(Dictionary), pointer :: isotope_dict - type(Dictionary), pointer :: ace_dict + ! These dictionaries provide a fast lookup mechanism + type(DictionaryII), pointer :: cell_dict + type(DictionaryII), pointer :: surface_dict + type(DictionaryII), pointer :: material_dict + type(DictionaryII), pointer :: isotope_dict + type(DictionaryCI), pointer :: xsdata_dict + type(DictionaryCI), pointer :: ace_dict ! Cross section arrays - type(AceContinuous), pointer :: xs_continuous(:) + type(AceContinuous), allocatable, target :: xs_continuous(:) type(AceThermal), allocatable, target :: xs_thermal(:) + integer :: n_continuous + integer :: n_thermal + + ! Current cell, surface, material + type(Cell), pointer :: cCell + type(Surface), pointer :: cSurface + type(Material), pointer :: cMaterial ! unionized energy grid - real(8), allocatable :: energy_grid(:) + integer :: n_grid ! number of points on unionized grid + real(8), allocatable :: e_grid(:) ! energies on unionized grid ! Histories/cycles/etc for both external source and criticality integer :: n_particles ! # of particles (per cycle for criticality) @@ -50,6 +60,7 @@ module global & MASS_NEUTRON = 1.0086649156, & ! mass of a neutron & MASS_PROTON = 1.00727646677, & ! mass of a proton & AMU = 1.66053873e-27, & ! 1 amu in kg + & N_AVOGADRO = 0.602214179, & ! Avogadro's number in 10^24/mol & K_BOLTZMANN = 8.617342e-5, & ! Boltzmann constant in eV/K & INFINITY = huge(0.0_8) ! positive infinity @@ -106,13 +117,13 @@ module global & ELECTRON_ = 3 integer :: verbosity = 5 - integer, parameter :: max_words = 100 + integer, parameter :: max_words = 500 integer, parameter :: max_line = 250 ! Versioning numbers integer, parameter :: VERSION_MAJOR = 0 integer, parameter :: VERSION_MINOR = 1 - integer, parameter :: VERSION_RELEASE = 2 + integer, parameter :: VERSION_RELEASE = 3 contains diff --git a/src/input_sample b/src/input_sample index 0988914339..4975c01d2e 100644 --- a/src/input_sample +++ b/src/input_sample @@ -1,15 +1,17 @@ # test input file -cell 1 1 -2 -cell 2 1 2 -1 +cell 100 40 -1 +cell 200 40 1 -2 -surface 1 sph 0 0 0 5 -surface 2 sph 0 0.0 0.0 3 +surface 1 sph 0 0 0 3 +surface 2 sph 0 0 0 5 -material 1 1001.03c 2.0 8016.03c 1.0 92235.12c 3.0 +material 40 -20.0 & + 3007.03c 1.0 & + 92238.03c 1.0 source box -3 -3 -3 3 3 3 xs_library endfb7 xs_data /opt/serpent/xsdata/endfb7 -criticality 1 1 10 +criticality 1 1 50 diff --git a/src/main.f90 b/src/main.f90 index 7370abdbc8..1a9fac72ab 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -4,15 +4,14 @@ program main use fileio, only: read_input, read_command_line, read_count, & & normalize_ao use output, only: title, echo_input, message, warning, error - use geometry, only: sense, cell_contains + use geometry, only: sense, cell_contains, neighbor_lists use mcnp_random, only: RN_init_problem, rang, RN_init_particle use source, only: init_source, get_source_particle use physics, only: transport - use data_structures, only: Dictionary, dict_create, dict_add_key, & - & dict_get_key - use cross_section, only: read_xsdata + use data_structures, only: dict_create, dict_add_key, dict_get_key + use cross_section, only: read_xsdata, material_total_xs use ace, only: read_xs - use energy_grid, only: unionized_grid + use energy_grid, only: unionized_grid, original_indices implicit none @@ -21,7 +20,7 @@ program main ! Print the OpenMC title and version/date/time information call title() - verbosity = 10 + verbosity = 9 ! Initialize random number generator call RN_init_problem( 3, 0_8, 0_8, 0_8, 0 ) @@ -37,20 +36,36 @@ program main ! pass to actually read values call read_count(path_input) call read_input(path_input) - call read_xsdata(path_xsdata) - call normalize_ao() - call read_xs() - call unionized_grid() - stop - + ! After reading input and basic geometry setup is complete, build + ! lists of neighboring cells for efficient tracking + call neighbor_lists() + + ! Read cross section summary file to determine what files contain + ! cross-sections + call read_xsdata(path_xsdata) + + ! With the AWRs from the xsdata, change all material specifications + ! so that they contain atom percents summing to 1 + call normalize_ao() + + ! Read ACE-format cross sections + call read_xs() + + ! Construct unionized energy grid from cross-sections + call unionized_grid() + call original_indices() + + ! calculate total material cross-sections for sampling path lenghts + call material_total_xs() + call echo_input() ! create source particles call init_source() ! start problem - surfaces(1)%bc = BC_VACUUM + surfaces(2)%bc = BC_VACUUM call run_problem() diff --git a/src/physics.f90 b/src/physics.f90 index eaf3db6005..2146e6f96d 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -5,6 +5,7 @@ module physics use types, only: Neutron use mcnp_random, only: rang use output, only: error, message + use search, only: binary_search implicit none @@ -15,7 +16,7 @@ contains ! geometry. !===================================================================== - subroutine transport( neut ) + subroutine transport(neut) type(Neutron), pointer, intent(inout) :: neut @@ -26,30 +27,56 @@ contains real(8) :: d_to_boundary real(8) :: d_to_collision real(8) :: distance + real(8) :: Sigma ! total cross-section + real(8) :: f ! interpolation factor + integer :: IE ! index on energy grid - if ( neut%cell == 0 ) then - call find_cell( neut ) + ! determine what cell the particle is in + if (neut%cell == 0) then + call find_cell(neut) + end if + if (verbosity >= 10) then + msg = "=== Particle " // trim(int_to_str(neut%uid)) // " ===" + call message(msg, 10) + + i = cells(neut%cell)%uid + msg = " Born in cell " // trim(int_to_str(i)) + call message(msg, 10) end if - do while ( neut%alive ) + ! sample energy from Watt fission energy spectrum for U-235 + neut%E = watt_spectrum(0.988_8, 2.249_8) - ! Determine distance neutron moves - call dist_to_boundary( neut, d_to_boundary, surf ) - d_to_collision = -log(rang()) / 1.0 - distance = min( d_to_boundary, d_to_collision ) + ! find energy index, interpolation factor + IE = binary_search(e_grid, n_grid, neut%E) + f = (neut%E - e_grid(IE))/(e_grid(IE+1) - e_grid(IE)) + + ! Determine material total cross-section + Sigma = f*cMaterial%total_xs(IE) + (1-f)*cMaterial%total_xs(IE+1) + neut%IE = IE + neut%interp = f + + do while (neut%alive) + + ! Find the distance to the nearest boundary + call dist_to_boundary(neut, d_to_boundary, surf) + + ! Sample a distance to collision + d_to_collision = -log(rang()) / 1.0 ! Sigma + distance = min(d_to_boundary, d_to_collision) ! Advance neutron neut%xyz = neut%xyz + distance*neut%uvw ! Add pathlength tallies - if ( d_to_collision > d_to_boundary ) then + if (d_to_collision > d_to_boundary) then neut%surface = surf neut%cell = 0 - call cross_boundary( neut ) + call cross_boundary(neut) else ! collision - call collision( neut ) + call collision(neut) end if end do @@ -60,20 +87,62 @@ contains ! COLLISION !===================================================================== - subroutine collision( neut ) + subroutine collision(neut) - type(Neutron), pointer, intent(inout) :: neut + type(Neutron), pointer :: neut + + type(AceContinuous), pointer :: table real(8) :: r1 real(8) :: phi ! azimuthal angle real(8) :: mu ! cosine of polar angle character(250) :: msg + integer :: cell_num + integer :: i,j + integer :: n_isotopes + integer :: IE + + real(8) :: f, Sigma + real(8) :: density, density_i + real(8) :: p + real(8), allocatable :: Sigma_t(:) ! tallies - + + density = cMaterial%atom_density + + ! calculate total cross-section for each nuclide at current energy + ! in order to create discrete pdf for sampling nuclide + n_isotopes = cMaterial%n_isotopes + allocate(Sigma_t(n_isotopes)) + do i = 1, n_isotopes + table => xs_continuous(cMaterial%table(i)) + density_i = cMaterial%atom_percent(i)*density + + ! search nuclide energy grid + IE = table%grid_index(neut%IE) + f = (neut%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE)) + + Sigma = density_i*(f*table%sigma_t(IE) + (1-f)*(table%sigma_t(IE+1))) + Sigma_t(i) = Sigma + end do + + ! normalize to create a discrete pdf + Sigma_t = Sigma_t / sum(Sigma_t) + + ! sample nuclide + r1 = rang() + p = 0.0_8 + do i = 1, n_isotopes + p = p + Sigma_t(i) + if (r1 < p) exit + end do + table => xs_continuous(cMaterial%table(i)) + ! print *, 'sampled nuclide ', i + ! select collision type r1 = rang() - if ( r1 <= 0.5 ) then + if (r1 <= 0.5) then ! scatter phi = 2.*pi*rang() mu = 2.*rang() - 1 @@ -82,12 +151,68 @@ contains neut%uvw(3) = sqrt(1. - mu**2) * sin(phi) else neut%alive = .false. - msg = "Particle " // trim(int_to_str(neut%uid)) // " was absorbed in cell " & - & // trim(int_to_str(neut%cell)) - call message( msg, 10 ) + if (verbosity >= 10) then + cell_num = cells(neut%cell)%uid + msg = " Absorbed in cell " // trim(int_to_str(cell_num)) + call message(msg, 10) + end if return end if end subroutine collision +!===================================================================== +! MAXWELL_SPECTRUM samples an energy from the Maxwell fission +! distribution based on a rejection sampling scheme. This is described +! in the MCNP manual volume I -- need to verify formula +!===================================================================== + + function maxwell_spectrum(T) result(E_out) + + real(8), intent(in) :: T ! tabulated function of incoming E + real(8) :: E_out ! sampled energy + + real(8) :: r1, r2, r3, r4 ! random numbers + real(8) :: d ! r1^2 + r2^2 + + r1 = rang() + do + r2 = rang() + d = r1*r1 + r2*r2 + if (d < 1) exit + r1 = r2 + end do + + r3 = rang() + r4 = rang() + E_out = -T*(r1**2 * log(r3) / d + log(r4)) + + end function maxwell_spectrum + +!===================================================================== +! WATT_SPECTRUM samples the outgoing energy from a Watt +! energy-dependent fission spectrum. Although fitted parameters exist +! for many nuclides, generally the continuous tabular distributions +! (LAW 4) should be used in lieu of the Watt spectrum +!===================================================================== + + function watt_spectrum(a, b) result(E_out) + + real(8), intent(in) :: a + real(8), intent(in) :: b + real(8) :: E_out + + real(8) :: g + real(8) :: r1, r2 + + g = sqrt((1 + a*b/8)**2 - 1) + (1 + a*b/8) + do + r1 = log(rang()) + r2 = log(rang()) + E_out = -a*g*r1 + if (((1 - g)*(1 - r1) - r2)**2 < b*E_out) exit + end do + + end function watt_spectrum + end module physics diff --git a/src/search.f90 b/src/search.f90 index 6763d593c3..48e9b4b8fa 100644 --- a/src/search.f90 +++ b/src/search.f90 @@ -1,5 +1,7 @@ module search + use output, only: error + contains !===================================================================== @@ -8,7 +10,7 @@ contains ! energy grid searching !===================================================================== - function binary_search( array, n, val ) result( index ) + function binary_search(array, n, val) result(index) real(8), intent(in) :: array(n) integer, intent(in) :: n @@ -18,12 +20,14 @@ contains integer :: L integer :: R real(8) :: testval + character(250) :: msg L = 1 R = n if (val < array(L) .or. val > array(R)) then - ! error + msg = "Value outside of array during binary search" + call error(msg) end if do while (R - L > 1) @@ -32,8 +36,8 @@ contains if (val > array(L) .and. val < array(L+1)) then index = L return - elseif (val > array(R+1) .and. val < array(R)) then - index = R + elseif (val > array(R-1) .and. val < array(R)) then + index = R-1 return end if @@ -41,9 +45,9 @@ contains index = L + (R - L)/2 testval = array(index) if (val > testval) then - L = index + 1 + L = index elseif (val < testval) then - R = index - 1 + R = index end if end do diff --git a/src/source.f90 b/src/source.f90 index 086e6e861e..1d73b45457 100644 --- a/src/source.f90 +++ b/src/source.f90 @@ -41,8 +41,8 @@ contains source_bank(i)%xyz = p_min + r*(p_max - p_min) ! angle - phi = 2.*PI*rang() - mu = 2.*rang() - 1 + phi = 2.0_8*PI*rang() + mu = 2.0_8*rang() - 1.0_8 source_bank(i)%uvw(1) = mu source_bank(i)%uvw(2) = sqrt(1. - mu**2) * cos(phi) source_bank(i)%uvw(3) = sqrt(1. - mu**2) * sin(phi) diff --git a/src/string.f90 b/src/string.f90 index 8dc2f0dc51..c867fe6d03 100644 --- a/src/string.f90 +++ b/src/string.f90 @@ -129,11 +129,11 @@ contains ! string = concatenated string !===================================================================== - subroutine concatenate( words, n_words, string ) + function concatenate(words, n_words) result(string) character(*), intent(in) :: words(n_words) integer, intent(in) :: n_words - character(250), intent(out) :: string + character(250) :: string integer :: i ! index @@ -143,7 +143,7 @@ contains string = trim(string) // ' ' // words(i) end do - end subroutine concatenate + end function concatenate !===================================================================== ! LOWER_CASE converts a string to all lower case characters diff --git a/src/types.f90 b/src/types.f90 index f3f9777af2..0c0d7cce61 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -33,13 +33,16 @@ module types !===================================================================== type Neutron - integer :: uid ! Unique ID - real(8) :: xyz(3) ! location - real(8) :: uvw(3) ! directional cosines - integer :: cell ! current cell - integer :: surface ! current surface - real(8) :: wgt ! particle weight - logical :: alive ! is particle alive? + integer :: uid ! Unique ID + real(8) :: xyz(3) ! location + real(8) :: uvw(3) ! directional cosines + real(8) :: E ! energy + integer :: IE ! index on energy grid + real(8) :: interp ! interpolation factor for energy grid + integer :: cell ! current cell + integer :: surface ! current surface + real(8) :: wgt ! particle weight + logical :: alive ! is particle alive? end type Neutron !===================================================================== @@ -74,10 +77,12 @@ module types type Material integer :: uid integer :: n_isotopes - character(10), allocatable :: names(:) - integer, allocatable :: isotopes(:) - integer, allocatable :: table(:) + character(10), allocatable :: names(:) ! isotope names + integer, allocatable :: isotopes(:) ! index in xsdata list + integer, allocatable :: table(:) ! index in xs array + real(8) :: atom_density ! total atom density in atom/b-cm real(8), allocatable :: atom_percent(:) + real(8), allocatable :: total_xs(:) ! macroscopic cross-section integer :: sab_table end type Material @@ -119,6 +124,8 @@ module types character(20) :: name real(8) :: awr real(8) :: temp + integer :: n_grid + integer, allocatable :: grid_index(:) real(8), allocatable :: energy(:) real(8), allocatable :: sigma_t(:) real(8), allocatable :: sigma_a(:) @@ -160,56 +167,6 @@ module types real(8), allocatable :: elastic_mu_out(:) end type AceThermal -!===================================================================== -! LISTDATA Data stored in a linked list. In this case, we store the -! (key,value) pair for a dictionary. Note that we need to store the -! key in addition to the value for collision resolution. -!===================================================================== - - ! Key length for dictionary - integer, parameter :: DICT_KEY_LENGTH = 20 - - type ListData - character(len=DICT_KEY_LENGTH) :: key - integer :: value - end type ListData - -!===================================================================== -! LINKEDLIST stores a simple linked list -!===================================================================== - - type LinkedList - type(LinkedList), pointer :: next - type(ListData) :: data - end type LinkedList - -!===================================================================== -! LINKEDLISTGRID stores a sorted list of energies for the unionized -! energy grid as a linked list -!===================================================================== - - type LinkedListGrid - type(LinkedListGrid), pointer :: next - real(8) :: energy - end type LinkedListGrid - -!===================================================================== -! HASHLIST - Since it's not possible to directly do an array of -! pointers, this derived type provides a pointer -!===================================================================== - - type HashList - type(LinkedList), pointer :: list - end type HashList - -!===================================================================== -! DICTIONARY provides a dictionary data structure of (key,value) pairs -!===================================================================== - - type Dictionary - type(HashList), pointer :: table(:) - end type Dictionary - !===================================================================== ! XSDATA contains data read in from a SERPENT xsdata file !===================================================================== @@ -226,4 +183,104 @@ module types character(100) :: path end type xsData +!===================================================================== +! KEYVALUECI stores the (key,value) pair for a dictionary where the +! key is a string and the value is an integer. Note that we need to +! store the key in addition to the value for collision resolution. +!===================================================================== + + ! Key length for dictionary + integer, parameter :: DICT_KEY_LENGTH = 20 + + type KeyValueCI + character(len=DICT_KEY_LENGTH) :: key + integer :: value + end type KeyValueCI + +!===================================================================== +! KEYVALUEII stores the (key,value) pair for a dictionary where the +! key is an integer and the value is an integer. Note that we need to +! store the key in addition to the value for collision resolution. +!===================================================================== + + type KeyValueII + integer :: key + integer :: value + end type KeyValueII + +!===================================================================== +! LISTKEYVALUECI stores a linked list of (key,value) pairs where the +! key is a character and the value is an integer +!===================================================================== + + type ListKeyValueCI + type(ListKeyValueCI), pointer :: next + type(KeyValueCI) :: data + end type ListKeyValueCI + +!===================================================================== +! LISTKEYVALUEII stores a linked list of (key,value) pairs where the +! key is a character and the value is an integer +!===================================================================== + + type ListKeyValueII + type(ListKeyValueII), pointer :: next + type(KeyValueII) :: data + end type ListKeyValueII + +!===================================================================== +! LISTREAL stores a linked list of real values. This is used for +! constructing a unionized energy grid. +!===================================================================== + + type ListReal + type(ListReal), pointer :: next + real(8) :: data + end type ListReal + +!===================================================================== +! LISTINT stores a linked list of integer values. +!===================================================================== + + type ListInt + type(ListInt), pointer :: next + integer :: data + end type ListInt + +!===================================================================== +! HASHLISTCI - Since it's not possible to directly do an array of +! pointers, this derived type provides a pointer +!===================================================================== + + type HashListCI + type(ListKeyValueCI), pointer :: list + end type HashListCI + +!===================================================================== +! HASHLISTII - Since it's not possible to directly do an array of +! pointers, this derived type provides a pointer +!===================================================================== + + type HashListII + type(ListKeyValueII), pointer :: list + end type HashListII + +!===================================================================== +! DICTIONARYCI provides a dictionary data structure of (key,value) +! pairs where the keys are strings and values are integers. +!===================================================================== + + type DictionaryCI + type(HashListCI), pointer :: table(:) + end type DictionaryCI + +!===================================================================== +! DICTIONARYII provides a dictionary data structure of (key,value) +! pairs where the keys and values are both integers. +!===================================================================== + + type DictionaryII + type(HashListII), pointer :: table(:) + end type DictionaryII + end module types