From bb61da69a2a7086f4a312561c1c6da25b76a8a4e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 27 Jan 2011 04:30:20 +0000 Subject: [PATCH] Added dictionary/linked list capability, other minor changes. --- ChangeLog | 10 + src/Makefile | 30 ++- src/data_structures.f90 | 495 ++++++++++++++++++++++++++++++++++++++++ src/geometry.f90 | 5 +- src/global.f90 | 5 +- src/main.f90 | 7 +- src/output.f90 | 12 +- src/physics.f90 | 39 +++- 8 files changed, 584 insertions(+), 19 deletions(-) create mode 100644 src/data_structures.f90 diff --git a/ChangeLog b/ChangeLog index a58ede0fd..fde6957f2 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,13 @@ +2011-01-26 Paul Romano + + * physics.f90: Added collision routine, still segfaults though. + * geometry.f90: Added messages for high verbosities + * global.f90: Added verbosity parameter, key length for + dictionaries + * output.f90: Implemented verbosity for messages. + * data_structures.f90: Added implementations of linked list and + dictionaries based on the flib open source package + 2011-01-23 Paul Romano * types.f90: Added cell, surface, and alive attributes to Neutron diff --git a/src/Makefile b/src/Makefile index 41d0ab1e0..b0c122fd8 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,8 +1,17 @@ program = openmc -src = main.f90 types.f90 global.f90 fileio.f90 output.f90 \ - string.f90 geometry.f90 mcnp_random.f90 source.f90 \ - physics.f90 +src = data_structures.f90 \ + fileio.f90 \ + geometry.f90 \ + global.f90 \ + main.f90 \ + mcnp_random.f90 \ + output.f90 \ + physics.f90 \ + source.f90 \ + string.f90 \ + types.f90 + objects = $(src:.f90=.o) #-------------------------------------------------------------------- @@ -30,13 +39,14 @@ neat: #-------------------------------------------------------------------- # Dependencies -types.o: -global.o: types.o -string.o: global.o output.o -output.o: global.o -source.o: global.o mcnp_random.o -physics.o: types.o global.o mcnp_random.o geometry.o output.o -geometry.o: types.o global.o output.o string.o +data_structures.o: global.o fileio.o: types.o global.o string.o output.o +geometry.o: types.o global.o output.o string.o +global.o: types.o main.o: global.o fileio.o output.o geometry.o mcnp_random.o \ source.o physics.o +output.o: global.o +physics.o: types.o global.o mcnp_random.o geometry.o output.o +source.o: global.o mcnp_random.o +string.o: global.o output.o +types.o: diff --git a/src/data_structures.f90 b/src/data_structures.f90 new file mode 100644 index 000000000..faca8b751 --- /dev/null +++ b/src/data_structures.f90 @@ -0,0 +1,495 @@ +module data_structures + +!===================================================================== +! DATA_STRUCTURES module +! +! This module implements a dictionary that has (key,value) pairs. This +! data structure is used to provide lookup features, e.g. cells and +! surfaces by name. +! +! The original version from the 'flibs' open source package +! implemented another derived type called DICT_DATA that has been +! replaced here by a simple integer in ListData. If used in the +! original form , the dictionary can store derived types (changes made +! to ListData, dict_create, dict_add_key, and dict_get_key). +!===================================================================== + + use global, only: DICT_KEY_LENGTH + + implicit none + + ! 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. + type ListData + character(len=DICT_KEY_LENGTH) :: key + integer :: value + end type ListData + + ! A simple linked list + type LinkedList + type(LinkedList), pointer :: next + type(ListData) :: data + end type LinkedList + + type HashList + type(LinkedList), pointer :: list + end type HashList + + ! A dictionary of (key,value) pairs + type Dictionary + private + type(HashList), pointer, dimension(:) :: table + end type Dictionary + + ! 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 + +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. +!===================================================================== + + subroutine list_create( list, data ) + + type(LinkedList), pointer :: list + type(ListData), intent(in) :: data + + allocate( list ) + list%next => null() + list%data = data + + end subroutine list_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 +!===================================================================== + + subroutine list_destroy( list ) + + type(LinkedList), pointer :: list + + type(LinkedList), pointer :: current + type(LinkedList), pointer :: next + + current => list + do while ( associated(current%next) ) + next => current%next + deallocate( current ) + current => next + enddo + + end subroutine list_destroy + +!===================================================================== +! LIST_COUNT count the number of items in the list +! Arguments: +! list Pointer to the list +!===================================================================== + + integer function list_count( list ) + + type(LinkedList), pointer :: list + + type(LinkedList), pointer :: current + type(LinkedList), pointer :: next + + if ( associated(list) ) then + list_count = 1 + current => list + do while ( associated(current%next) ) + current => current%next + list_count = list_count + 1 + enddo + else + list_count = 0 + endif + + end function list_count + +!===================================================================== +! LIST_NEXT returns the next element (if any) +! Arguments: +! elem Element in the linked list +! Result: +!===================================================================== + + function list_next( elem ) result(next) + + type(LinkedList), pointer :: elem + type(LinkedList), pointer :: next + + next => elem%next + + end function list_next + +!===================================================================== +! LIST_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 ) + + type(LinkedList), pointer :: elem + type(ListData), intent(in) :: data + + type(LinkedList), pointer :: next + + allocate(next) + + next%next => elem%next + elem%next => next + next%data = data + + end subroutine list_insert + +!===================================================================== +! LIST_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 ) + + type(LinkedList), pointer :: list + type(ListData), intent(in) :: data + + type(LinkedList), pointer :: elem + + allocate(elem) + elem%data = data + + elem%next => list + list => elem + + end subroutine list_insert_head + +!===================================================================== +! LIST_DELETE_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 ) + + type(LinkedList), pointer :: list + type(LinkedList), pointer :: elem + + type(LinkedList), pointer :: current + type(LinkedList), pointer :: prev + + 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_delete_element + +!===================================================================== +! LIST_GET_DATA gets the data stored with a list element +! Arguments: +! elem Element in the linked list +!===================================================================== + + function list_get_data( elem ) result(data) + + type(LinkedList), pointer :: elem + type(ListData) :: data + + data = elem%data + + end function list_get_data + +!===================================================================== +! LIST_PUT_DATA stores new data with a list element +! Arguments: +! elem Element in the linked list +! data The data to be stored +!===================================================================== + + subroutine list_put_data( elem, data ) + + type(LinkedList), pointer :: elem + type(ListData), intent(in) :: data + + elem%data = data + + end subroutine list_put_data + +!===================================================================== +! DICT_CREATE creates and initializes a dictionary +! 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. +!===================================================================== + + subroutine dict_create( dict, key, value ) + + type(Dictionary), pointer :: dict + character(len=*), intent(in) :: key + integer, intent(in) :: value + + type(ListData) :: data + integer :: i + integer :: hash + + 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 + +!===================================================================== +! 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 +!===================================================================== + + subroutine dict_destroy( dict ) + + type(Dictionary), pointer :: dict + + integer :: i + + do i = 1,size(dict%table) + if ( associated( dict%table(i)%list ) ) then + call list_destroy( dict%table(i)%list ) + endif + enddo + deallocate( dict%table ) + deallocate( dict ) + + end subroutine dict_destroy + +!===================================================================== +! DICT_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_add_key( dict, key, value ) + + type(Dictionary), pointer :: dict + character(len=*), intent(in) :: key + integer, intent(in) :: value + + type(ListData) :: data + type(LinkedList), 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( 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 ) + endif + endif + + end subroutine dict_add_key + +!===================================================================== +! DICT_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 ) + + type(Dictionary), pointer :: dict + character(len=*), intent(in) :: key + + type(LinkedList), pointer :: elem + integer :: hash + + elem => dict_get_elem( dict, key ) + + if ( associated(elem) ) then + hash = dict_hashkey( trim(key) ) + call list_delete_element( dict%table(hash)%list, elem ) + endif + + end subroutine dict_delete_key + +!===================================================================== +! DICT_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) + + type(Dictionary), pointer :: dict + character(len=*), intent(in) :: key + integer :: value + + type(ListData) :: data + type(LinkedList), pointer :: elem + + elem => dict_get_elem( dict, key ) + + if ( associated(elem) ) then + value = elem%data%value + else + value = DICT_NULL + endif + + end function dict_get_key + +!===================================================================== +! DICT_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) + + type(Dictionary), pointer :: dict + character(len=*), intent(in) :: key + logical :: has + + type(LinkedList), pointer :: elem + + elem => dict_get_elem( dict, key ) + + has = associated(elem) + + end function dict_has_key + +!===================================================================== +! DICT_GET_ELEM finds the element with a particular key +! Arguments: +! dict Pointer to the dictionary +! key Key to be sought +!===================================================================== + + function dict_get_elem( dict, key ) result(elem) + + type(Dictionary), pointer :: dict + character(len=*), intent(in) :: key + + type(LinkedList), 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 + exit + else + elem => list_next( elem ) + endif + enddo + + end function dict_get_elem + +!===================================================================== +! DICT_HASHKEY determine the hash value from the string +! Arguments: +! key String to be examined +!===================================================================== + + integer function dict_hashkey( key ) + + character(len=*), intent(in) :: key + + integer :: hash + integer :: i + + dict_hashkey = 0 + + do i = 1,len(key) + dict_hashkey = multiplier * dict_hashkey + 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 ) + + end function dict_hashkey + +end module data_structures diff --git a/src/geometry.f90 b/src/geometry.f90 index 34dee79cb..9f941cffc 100644 --- a/src/geometry.f90 +++ b/src/geometry.f90 @@ -2,7 +2,7 @@ module geometry use global use types, only: Cell, Surface - use output, only: error + use output, only: error, message implicit none @@ -131,6 +131,9 @@ contains ! check for leakage 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 ) return end if diff --git a/src/global.f90 b/src/global.f90 index 5d3a720ac..03c865bdc 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -83,7 +83,7 @@ module global character(32) :: inputfile - integer, parameter :: verbosity = 5 + integer :: verbosity = 5 integer, parameter :: max_words = 100 ! Versioning numbers @@ -91,6 +91,9 @@ module global integer, parameter :: VERSION_MINOR = 1 integer, parameter :: VERSION_RELEASE = 1 + ! Key length for dictionary + integer, parameter :: DICT_KEY_LENGTH = 20 + contains !===================================================================== diff --git a/src/main.f90 b/src/main.f90 index c837d8438..d301fbc65 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -7,15 +7,14 @@ program main 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 implicit none character(16) :: filename character(250) :: msg - real(8) :: point(3) - integer :: s(4) - ! Print the OpenMC title and version/date/time information call title() @@ -39,6 +38,8 @@ program main call init_source() ! start problem + verbosity = 10 + surfaces(1)%bc = BC_VACUUM call run_problem() diff --git a/src/output.f90 b/src/output.f90 index d84da2697..d08b7fd10 100644 --- a/src/output.f90 +++ b/src/output.f90 @@ -102,8 +102,18 @@ module output character(*), intent(in) :: msg integer, intent(in) :: level + integer :: ou + integer :: n_lines + integer :: i + if ( level <= verbosity ) then - write (OUTPUT_UNIT,*) trim(msg) + ou = OUTPUT_UNIT + + n_lines = (len_trim(msg)-1)/79 + 1 + do i = 1, n_lines + write(ou, fmt='(1X,A79)') msg(79*(i-1)+1:79*i) + end do + end if end subroutine message diff --git a/src/physics.f90 b/src/physics.f90 index c9b044aef..eaf3db600 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -4,7 +4,7 @@ module physics use geometry, only: find_cell, dist_to_boundary, cross_boundary use types, only: Neutron use mcnp_random, only: rang - use output, only: error + use output, only: error, message implicit none @@ -49,12 +49,45 @@ contains call cross_boundary( neut ) else ! collision - msg = "Collision not implemented yet!" - call error( msg ) + call collision( neut ) end if end do end subroutine transport +!===================================================================== +! COLLISION +!===================================================================== + + subroutine collision( neut ) + + type(Neutron), pointer, intent(inout) :: neut + + real(8) :: r1 + real(8) :: phi ! azimuthal angle + real(8) :: mu ! cosine of polar angle + character(250) :: msg + + ! tallies + + ! select collision type + r1 = rang() + if ( r1 <= 0.5 ) then + ! scatter + phi = 2.*pi*rang() + mu = 2.*rang() - 1 + neut%uvw(1) = mu + neut%uvw(2) = sqrt(1. - mu**2) * cos(phi) + 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 ) + return + end if + + end subroutine collision + end module physics