crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Read ACE continuous-energy tables, energy grid unionization, reading materials, other changes.

Browse source
This commit is contained in:
Paul Romano 2011-02-03 05:47:35 +00:00
parent bb61da69a2
commit e64b4f53b9
15 changed files with 1349 additions and 169 deletions
Download patch file Download diff file Expand all files Collapse all files

24
ChangeLog
View file

@ -1,3 +1,27 @@
2011-02-03 Paul Romano <romano7@mit.edu>
* ace.f90: Added initial routines for reading continuous-energy
ACE cross section tables. No secondary energy/angle distributions
or nu data is read.
* cross_section.f90: Added read_xsdata routine which reads
information from an xsdata file from SERPENT.
* data_structures.f90: Moved derived types to the types module and
change dict_create so that one doesn't need to add a (key,value)
pair on initialization.
* energy_grid.f90: First routines for energy grid unionization
based on insertion sorts.
* fileio.f90: Added routines for reading integer and real data,
particularly from cross-section libraries. Made an interface so
that a single routine can be called regardless of type. Ability to
read material cards and added subroutine for normalization of
atom/weight percents.
* Makefile: Added 'unittest' target
* output.f90: Fixed bug in get_today subroutine whereby the day
was displayed wrong.
* search.f90: Added binary search which will eventually be used
for energy grid searches
* unittest.f90: Added this program for unit testing.
2011-01-26 Paul Romano <romano7@mit.edu>
* physics.f90: Added collision routine, still segfaults though.

30
src/Makefile
View file

@ -1,18 +1,24 @@
program = openmc
src = data_structures.f90 \
main = main.f90
test = unittest.f90
modules = ace.f90 \
cross_section.f90 \
data_structures.f90 \
energy_grid.f90 \
fileio.f90 \
geometry.f90 \
global.f90 \
main.f90 \
mcnp_random.f90 \
output.f90 \
physics.f90 \
search.f90 \
source.f90 \
string.f90 \
types.f90
types.f90 \
objects = $(src:.f90=.o)
main_objects = $(modules:.f90=.o) $(main:.f90=.o)
test_objects = $(modules:.f90=.o) $(test:.f90=.o)
#--------------------------------------------------------------------
# Compiler Options
@ -22,12 +28,14 @@ F90FLAGS = -g
#--------------------------------------------------------------------
# Targets
all: $(program)
$(program): $(objects)
$(F90) $(F90FLAGS) $(objects) -o $@
$(program): $(main_objects)
$(F90) $(F90FLAGS) $(main_objects) -o $@
clean:
@rm -f *.o *.mod $(program)
neat:
@rm -f *.o *.mod
unittest: $(test_objects)
$(F90) $(F90FLAGS) $(test_objects) -o $@
#--------------------------------------------------------------------
# Rules & misc
@ -39,14 +47,20 @@ neat:
#--------------------------------------------------------------------
# Dependencies
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
fileio.o: types.o global.o string.o output.o
energy_grid.o: global.o
fileio.o: types.o global.o string.o output.o data_structures.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
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:
source.o: global.o mcnp_random.o
string.o: global.o output.o
types.o:
unittest.o: global.o energy_grid.o

328
src/ace.f90 Normal file
View file

@ -0,0 +1,328 @@
module ace
use global
use output, only: error, message
use string, only: lower_case
use fileio, only: read_line, read_data, skip_lines
use string, only: split_string, str_to_real
use data_structures, only: dict_create, dict_add_key, dict_has_key, &
& dict_get_key
integer :: NXS(16)
integer :: JXS(32)
real(8), allocatable :: XSS(:)
integer :: XSS_index
private :: NXS
private :: JXS
private :: XSS
contains
!=====================================================================
! READ_XS reads all the cross sections for the problem and stores them
! in xs_continuous and xs_thermal arrays
!=====================================================================
subroutine read_xs()
type(Material), pointer :: mat => null()
type(xsData), pointer :: iso => null()
type(AceContinuous), pointer :: ace_cont => null()
type(AceThermal), pointer :: ace_thermal => null()
integer :: i, j
integer :: index
character(10) :: table
character(250) :: msg
integer :: n
integer :: n_continuous
integer :: n_thermal
integer :: index_continuous
integer :: index_thermal
! 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))
case ('c')
if (.not. dict_has_key(ace_dict, table)) then
index_continuous = index_continuous + 1
call dict_add_key(ace_dict, table, index_continuous)
call read_ACE_continuous(index_continuous, index)
mat%table(j) = index_continuous
end if
case ('t')
index_thermal = index_thermal + 1
n_thermal = n_thermal + 1
case default
msg = "Unknown cross section table type: " // table
call error(msg)
end select
end do
end do
end subroutine read_xs
!=====================================================================
! READ_ACE_CONTINUOUS reads in a single ACE continuous-energy cross
! section table
!=====================================================================
subroutine read_ACE_continuous(index_table, index)
integer, intent(in) :: index_table
integer, intent(in) :: index
type(AceContinuous), pointer :: table => null()
integer :: i
integer :: in = 7
integer :: ioError
integer :: words_per_line
integer :: lines
integer :: n
logical :: file_exists
logical :: found_xs
character(7) :: readable
character(250) :: msg, line
character(32) :: words(max_words)
character(100) :: filename
character(10) :: tablename
real(8) :: kT
filename = xsdatas(index)%path
tablename = xsdatas(index)%id
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
msg = "ACE library '" // trim(filename) // "' does not exist!"
call error( msg )
elseif ( readable(1:3) /= 'YES' ) then
msg = "ACE library '" // trim(filename) // "' is not readable! &
&Change file permissions with chmod command."
call error( msg )
end if
! display message
msg = "Loading ACE cross section table: " // tablename
call message(msg, 6)
! open file
open(file=filename, unit=in, status='old', action='read')
found_xs = .false.
do while (.not. found_xs)
call read_line(in, line, ioError)
if (ioError < 0) then
msg = "Could not find ACE table " // tablename // "."
call error(msg)
end if
call split_string(line, words, n)
if (trim(words(1)) == trim(tablename)) then
found_xs = .true.
table%name = words(1)
table%awr = str_to_real(words(2))
kT = str_to_real(words(3))
table%temp = kT / K_BOLTZMANN
end if
! Skip 5 lines
call skip_lines(in, 5, ioError)
! Read NXS data
lines = 2
words_per_line = 8
call read_data(in, NXS, 16, lines, words_per_line)
! Read JXS data
lines = 4
call read_data(in, JXS, 32, lines, words_per_line)
! Calculate how many data points and lines in the XSS array
n = NXS(1)
lines = (n + 3)/4
if (found_xs) then
! allocate storage for XSS array
allocate(XSS(n))
! Read XSS
words_per_line = 4
call read_data(in, XSS, n, lines, words_per_line)
else
call skip_lines(in, lines, ioError)
end if
end do
call parse_ESZ(table)
call parse_MTR(table)
! Free memory from XSS array
if(allocated(XSS)) deallocate(XSS)
if(associated(table)) nullify(table)
close(unit=in)
end subroutine read_ACE_continuous
!=====================================================================
! PARSE_ESZ - reads through the ESZ block. This block contains the
! energy grid, total xs, absorption xs, elastic scattering xs, and
! heating numbers.
!=====================================================================
subroutine parse_ESZ(table)
type(AceContinuous), pointer, intent(inout) :: table
integer :: NE
! determine number of energy points
NE = NXS(3)
! allocate storage for arrays
allocate(table%energy(NE))
allocate(table%sigma_t(NE))
allocate(table%sigma_a(NE))
allocate(table%sigma_el(NE))
allocate(table%heating(NE))
! read data from XSS
XSS_index = 1
table%energy = get_real(NE)
table%sigma_t = get_real(NE)
table%sigma_a = get_real(NE)
table%sigma_el = get_real(NE)
table%heating = get_real(NE)
end subroutine parse_ESZ
!=====================================================================
! PARSE_MTR - Get the list of reaction MTs for this cross-section
! table. The MT values are somewhat arbitrary. Also read in Q-values,
! neutron multiplicities, and cross-sections.
! =====================================================================
subroutine parse_MTR(table)
type(AceContinuous), pointer :: table
type(AceReaction), pointer :: rxn
integer :: LMT ! index of MT list in XSS
integer :: NMT ! Number of reactions
integer :: JXS4 ! index of Q values in XSS
integer :: JXS5 ! index of neutron multiplicities in XSS
integer :: JXS7 ! index of reactions cross-sections in XSS
integer :: LXS !
integer :: LOCA !
integer :: NE ! number of energies for reaction
LMT = JXS(3)
JXS4 = JXS(4)
JXS5 = JXS(5)
LXS = JXS(6)
JXS7 = JXS(7)
NMT = NXS(4)
! allocate array of reactions
allocate(table%reactions(NMT))
do i = 1, NMT
rxn => table%reactions(i)
! read MT number, Q-value, and neutrons produced
rxn%MT = XSS(LMT+i-1)
rxn%Q_value = XSS(JXS4+i-1)
rxn%TY = XSS(JXS5+i-1)
! read cross section values
LOCA = XSS(LXS+i-1)
rxn%energy_index = XSS(JXS7 + LOCA - 1)
NE = XSS(JXS7 + LOCA)
allocate(rxn%sigma(NE))
XSS_index = JXS7 + LOCA + 1
rxn%sigma = get_real(NE)
end do
end subroutine parse_MTR
!=====================================================================
! GET_INT returns an array of integers read from the current position
! in the XSS array
!=====================================================================
function get_int(n_values) result(array)
integer, intent(in) :: n_values
integer :: array(n_values)
array = int(XSS(XSS_index:XSS_index + n_values - 1))
XSS_index = XSS_index + n_values
end function get_int
!=====================================================================
! GET_REAL returns an array of real(8)s read from the current position
! in the XSS array
!=====================================================================
function get_real(n_values) result(array)
integer, intent(in) :: n_values
real(8) :: array(n_values)
array = XSS(XSS_index:XSS_index + n_values - 1)
XSS_index = XSS_index + n_values
end function get_real
!=====================================================================
! READ_ACE_THERMAL reads in a single ACE S(a,b) thermal scattering
! cross section table
!=====================================================================
subroutine read_ACE_thermal( filename )
character(*), intent(in) :: filename
end subroutine read_ACE_thermal
end module ace

113
src/cross_section.f90 Normal file
View file

@ -0,0 +1,113 @@
module cross_section
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 types, only: xsData
implicit none
contains
!=====================================================================
! READ_XSDATA reads the data in a SERPENT xsdata file and builds a
! dictionary to find cross-section information later on.
!=====================================================================
subroutine read_xsdata( path )
character(*), intent(in) :: path
type(xsData), pointer :: iso => null()
character(250) :: line, msg
character(100) :: words(max_words)
character(100) :: filename
integer :: i, n
integer :: in = 7
logical :: file_exists
character(7) :: readable
integer :: count
integer :: index
integer :: ioError
! Construct filename
n = len(path)
if (path(n:n) == '/') then
filename = trim(path) // 'xsdata'
else
filename = trim(path) // '/xsdata'
end if
! Check if xsdata exists and is readable
inquire( FILE=filename, EXIST=file_exists, READ=readable )
if ( .not. file_exists ) then
msg = "Cross section summary '" // trim(filename) // "' does not exist!"
call error( msg )
elseif ( readable(1:3) /= 'YES' ) 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')
! determine how many lines
count = 0
do
read(unit=in, fmt='(A250)', iostat=ioError) line
if ( ioError < 0 ) then
! reached end of file
exit
elseif ( ioError > 0 ) then
msg = "Unknown error while reading file: " // filename
close(unit=in)
call error( msg )
end if
count = count + 1
end do
allocate(xsdatas(count))
! create dictionary
call dict_create(xsdata_dict)
! read actual lines
index = 0
rewind(in)
do
read(unit=in, fmt='(A250)', iostat=ioError) line
if ( ioError < 0 ) exit
index = index + 1
call split_string(line, words, n)
if ( n == 0 ) cycle ! skip blank line
! Check to make sure there are enough arguments
if ( n < 9 ) then
msg = "Not enough arguments on xsdata line: " // line
close(unit=in)
call error( msg )
end if
iso => xsdatas(index)
! store data
iso%alias = words(1)
iso%id = words(2)
iso%type = str_to_int(words(3))
iso%zaid = str_to_int(words(4))
iso%isomeric = str_to_int(words(5))
iso%awr = str_to_real(words(6))
iso%temp = str_to_real(words(7))
iso%binary = str_to_int(words(8))
iso%path = words(9)
! create dictionary entry
call dict_add_key(xsdata_dict, iso%alias, index)
end do
close(unit=in)
end subroutine read_xsdata
end module cross_section

46
src/data_structures.f90
View file

@ -14,34 +14,10 @@ module data_structures
! to ListData, dict_create, dict_add_key, and dict_get_key).
!=====================================================================
use global, only: DICT_KEY_LENGTH
use types, only: ListData, LinkedList, HashList, Dictionary
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
@ -281,15 +257,15 @@ contains
! destroy up an old list.
!=====================================================================
subroutine dict_create( dict, key, value )
subroutine dict_create( dict ) !, key, value )
type(Dictionary), pointer :: dict
character(len=*), intent(in) :: key
integer, intent(in) :: value
!!$ character(len=*), intent(in) :: key
!!$ integer, intent(in) :: value
type(ListData) :: data
!!$ type(ListData) :: data
integer :: i
integer :: hash
!!$ integer :: hash
allocate( dict )
allocate( dict%table(hash_size) )
@ -298,11 +274,11 @@ contains
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 )
!!$ data%key = key
!!$ data%value = value
!!$
!!$ hash = dict_hashkey( trim(key ) )
!!$ call list_create( dict%table(hash)%list, data )
end subroutine dict_create

167
src/energy_grid.f90 Normal file
View file

@ -0,0 +1,167 @@
module energy_grid
use global
use output, only: message
contains
!=====================================================================
! UNIONIZED_GRID creates a single unionized energy grid combined from
! each nuclide of each material. Right now, the grid for each nuclide
! is added into a linked list one at a time with an effective
! insertion sort. Could be done with a hash for all energy points and
! then a quicksort at the end (what hash function to use?)
!=====================================================================
subroutine unionized_grid()
type(LinkedListGrid), pointer :: list => null()
type(LinkedListGrid), pointer :: current => null()
type(Material), pointer :: mat => null()
type(AceContinuous), pointer :: table => null()
type(AceReaction), pointer :: rxn => null()
integer :: i, j
integer :: n
character(100) :: msg
msg = "Creating unionized energy grid"
call message(msg, 5)
! loop over all materials
do i = 1, n_materials
mat => materials(i)
! loop over all isotopes
do j = 1, mat%n_isotopes
table => xs_continuous(mat%table(j))
! loop over energy points
n = size(table%energy)
call add_grid_points(list, table%energy, n)
end do
end do
end subroutine unionized_grid
!=====================================================================
! ADD_GRID_POINTS adds energy points from the 'energy' array into a
! linked list of points already stored from previous arrays.
!=====================================================================
subroutine add_grid_points(list, energy, n_energy)
type(LinkedListGrid), 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()
integer :: index
real(8) :: E
index = 1
! if the original list is empty, we need to allocate the first
! element and store first energy point
if (grid_count(list) == 0) then
allocate(list)
current => list
do index = 1, n_energy
current%energy = energy(index)
if (index == n_energy) then
current%next => null()
return
end if
allocate(current%next)
current => current%next
end do
end if
current => list
head => list
E = energy(index)
do while (index <= n_energy)
! If we've reached the end of the grid energy list, add the
! remaining energy points to the end
if (.not. associated(current)) then
! finish remaining energies
do while (index <= n_energy)
allocate(previous%next)
current => previous%next
current%energy = energy(index)
previous => current
index = index + 1
end do
current%next => null()
exit
end if
if (E < current%energy) then
! create new element and insert it in energy grid list
allocate(tmp)
tmp%energy = E
tmp%next => current
if (associated(previous)) then
previous%next => tmp
previous => tmp
else
previous => tmp
head => previous
end if
nullify(tmp)
! advance index
index = index + 1
E = energy(index)
elseif (E == current%energy) then
! found the exact same energy, no need to store duplicates
! so just skip and move to next index
index = index + 1
E = energy(index)
else
previous => current
current => current%next
end if
end do
! It's possible that an element was inserted at the front of the
! list, so we need to move the list pointer back to the start of
! the list
list => head
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.
!=====================================================================
function grid_count(list) result(count)
type(LinkedListGrid), pointer :: list
integer :: count
type(LinkedListGrid), pointer :: current
! 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
end function grid_count
end module energy_grid

385
src/fileio.f90
View file

@ -4,9 +4,21 @@ module fileio
use types, only: Cell, Surface, ExtSource
use string, only: split_string_wl, lower_case, str_to_real, split_string
use output, only: message, warning, error
use data_structures, only: Dictionary, dict_create, dict_add_key, dict_get_key
implicit none
!=====================================================================
! READ_DATA interface allows data to be read with one function
! regardless of whether it is integer or real data. E.g. NXS and JXS
! can be read with the integer version and XSS can be read with the
! real version
!=====================================================================
interface read_data
module procedure read_data_int, read_data_real
end interface
contains
!=====================================================================
@ -22,21 +34,21 @@ contains
argc = COMMAND_ARGUMENT_COUNT()
if ( argc > 0 ) then
call GET_COMMAND_ARGUMENT(1,inputfile)
call GET_COMMAND_ARGUMENT(1, path_input)
else
msg = "No input file specified!"
call error( msg )
call error(msg)
end if
! Check if input file exists and is readable
inquire( FILE=inputfile, EXIST=file_exists, READ=readable )
if ( .not. file_exists ) then
msg = "Input file '" // trim(inputfile) // "' does not exist!"
call error( msg )
elseif ( readable(1:3) /= 'YES' ) then
msg = "Input file '" // trim(inputfile) // "' is not readable! &
inquire(FILE=path_input, EXIST=file_exists, READ=readable)
if (.not. file_exists) then
msg = "Input file '" // trim(path_input) // "' does not exist!"
call error(msg)
elseif (readable(1:3) /= 'YES') then
msg = "Input file '" // trim(path_input) // "' is not readable! &
&Change file permissions with chmod command."
call error( msg )
call error(msg)
end if
end subroutine read_command_line
@ -65,40 +77,49 @@ contains
do
read(unit=in, fmt='(A250)', iostat=readError) line
if ( readError /= 0 ) exit
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)
if ( n == 0 ) cycle
if (n == 0) cycle
select case ( trim(words(1)) )
case ( 'cell' )
select case (trim(words(1)))
case ('cell')
n_cells = n_cells + 1
case ( 'surface' )
case ('surface')
n_surfaces = n_surfaces + 1
case ( 'material' )
case ('material')
n_materials = n_materials + 1
end select
end do
! Check to make sure there are cells, surface, and materials
! defined for the problem
if ( n_cells == 0 ) then
if (n_cells == 0) then
msg = "No cells specified!"
call error( msg )
elseif ( n_surfaces == 0 ) then
close(unit=in)
call error(msg)
elseif (n_surfaces == 0) then
msg = "No surfaces specified!"
call error( msg )
elseif ( n_materials == 0 ) then
close(unit=in)
call error(msg)
elseif (n_materials == 0) then
msg = "No materials specified!"
call error( msg )
close(unit=in)
call error(msg)
end if
close(unit=in)
! Allocate arrays for cells, surfaces, etc.
allocate( cells(n_cells) )
allocate( surfaces(n_surfaces) )
allocate( materials(n_materials) )
allocate(cells(n_cells))
allocate(surfaces(n_surfaces))
allocate(materials(n_materials))
! Set up dictionaries
call dict_create(cell_dict)
call dict_create(surface_dict)
end subroutine read_count
!=====================================================================
@ -134,26 +155,31 @@ contains
call split_string_wl(line, words, n)
! skip comments
if ( n==0 .or. words(1)(1:1) == '#' ) cycle
if (n==0 .or. words(1)(1:1) == '#') cycle
select case ( trim(words(1)) )
select case (trim(words(1)))
case ('cell')
! Read data from cell entry
index_cell = index_cell + 1
call read_cell( index_cell, words, n )
call read_cell(index_cell, words, n)
case ('surface')
! Read data
index_surface = index_surface + 1
call read_surface( index_surface, words, n )
call read_surface(index_surface, words, n)
case ( 'material' )
case ('material')
index_material = index_material + 1
call read_material(index_material, words, n)
case ( 'source' )
call read_source( words, n )
case ('source')
call read_source(words, n)
case ( 'criticality' )
call read_criticality( words, n )
case ('criticality')
call read_criticality(words, n)
case ('xs_data')
path_xsdata = words(2)
case default
! do nothing
@ -163,13 +189,15 @@ contains
close(unit=in)
! call adjust_cell_index
end subroutine read_input
!=====================================================================
! READ_CELL parses the data on a cell card.
!=====================================================================
subroutine read_cell( index, words, n_words )
subroutine read_cell(index, words, n_words)
integer, intent(in) :: index
character(*), intent(in) :: words(max_words)
@ -185,25 +213,26 @@ contains
! Read cell identifier
read(words(2), fmt='(I8)', iostat=readError) this_cell%uid
if ( readError > 0 ) then
if (readError > 0) then
msg = "Invalid cell name: " // words(2)
call error( msg )
call error(msg)
end if
call dict_add_key(cell_dict, words(2), index)
! Read cell material
read(words(3), fmt='(I8)', iostat=readError) this_cell%material
! Read list of surfaces
allocate( this_cell%boundary_list(n_words-3) )
allocate(this_cell%boundary_list(n_words-3))
do i = 1, n_words-3
word = words(i+3)
if ( word(1:1) == '(' ) then
if (word(1:1) == '(') then
this_cell%boundary_list(i) = OP_LEFT_PAREN
elseif ( word(1:1) == ')' ) then
elseif (word(1:1) == ')') then
this_cell%boundary_list(i) = OP_RIGHT_PAREN
elseif ( word(1:1) == ':' ) then
elseif (word(1:1) == ':') then
this_cell%boundary_list(i) = OP_UNION
elseif ( word(1:1) == '#' ) then
elseif (word(1:1) == '#') then
this_cell%boundary_list(i) = OP_DIFFERENCE
else
read(word, fmt='(I8)', iostat=readError) this_cell%boundary_list(i)
@ -216,7 +245,7 @@ contains
! READ_SURFACE parses the data on a surface card.
!=====================================================================
subroutine read_surface( index, words, n_words )
subroutine read_surface(index, words, n_words)
integer, intent(in) :: index
character(*), intent(in) :: words(max_words)
@ -233,67 +262,67 @@ contains
! Read surface identifier
read(words(2), fmt='(I8)', iostat=readError) this_surface%uid
if ( readError > 0 ) then
if (readError > 0) then
msg = "Invalid surface name: " // words(2)
call error( msg )
call error(msg)
end if
! Read surface type
word = words(3)
call lower_case(word)
select case ( trim(word) )
case ( 'px' )
select case (trim(word))
case ('px')
this_surface%type = SURF_PX
coeffs_reqd = 1
case ( 'py' )
case ('py')
this_surface%type = SURF_PY
coeffs_reqd = 1
case ( 'pz' )
case ('pz')
this_surface%type = SURF_PZ
coeffs_reqd = 1
case ( 'plane' )
case ('plane')
this_surface%type = SURF_PLANE
coeffs_reqd = 4
case ( 'cylx' )
case ('cylx')
this_surface%type = SURF_CYL_X
coeffs_reqd = 3
case ( 'cyly' )
case ('cyly')
this_surface%type = SURF_CYL_Y
coeffs_reqd = 3
case ( 'cylz' )
case ('cylz')
this_surface%type = SURF_CYL_Z
coeffs_reqd = 3
case ( 'sph' )
case ('sph')
this_surface%type = SURF_SPHERE
coeffs_reqd = 4
case ( 'boxx' )
case ('boxx')
this_surface%type = SURF_BOX_X
coeffs_reqd = 4
case ( 'boxy' )
case ('boxy')
this_surface%type = SURF_BOX_Y
coeffs_reqd = 4
case ( 'boxz' )
case ('boxz')
this_surface%type = SURF_BOX_Z
coeffs_reqd = 4
case ( 'box' )
case ('box')
this_surface%type = SURF_BOX
coeffs_reqd = 6
case ( 'gq' )
case ('gq')
this_surface%type = SURF_GQ
coeffs_reqd = 10
case default
msg = "Invalid surface type: " // words(3)
call error( msg )
call error(msg)
end select
! Make sure there are enough coefficients for surface type
if ( n_words-3 < coeffs_reqd ) then
if (n_words-3 < coeffs_reqd) then
msg = "Not enough coefficients for surface: " // words(2)
call error( msg )
call error(msg)
end if
! Read list of surfaces
allocate( this_surface%coeffs(n_words-3) )
allocate(this_surface%coeffs(n_words-3))
do i = 1, n_words-3
this_surface%coeffs(i) = str_to_real(words(i+3))
end do
@ -304,7 +333,7 @@ contains
! READ_SOURCE parses the data on a source card.
!=====================================================================
subroutine read_source( words, n_words )
subroutine read_source(words, n_words)
character(*), intent(in) :: words(max_words)
integer, intent(in) :: n_words
@ -318,29 +347,135 @@ contains
! Read source type
word = words(2)
call lower_case(word)
select case ( trim(word) )
case ( 'box' )
select case (trim(word))
case ('box')
external_source%type = SRC_BOX
values_reqd = 6
case default
msg = "Invalid source type: " // words(2)
call error( msg )
call error(msg)
end select
! Make sure there are enough values for this source type
if ( n_words-2 < values_reqd ) then
if (n_words-2 < values_reqd) then
msg = "Not enough values for source of type: " // words(2)
call error( msg )
call error(msg)
end if
! Read list of surfaces
allocate( external_source%values(n_words-2) )
allocate(external_source%values(n_words-2))
do i = 1, n_words-2
external_source%values(i) = str_to_real(words(i+2))
end do
end subroutine read_source
!=====================================================================
! READ_MATERIAL parses a material card, creating the appropriate
! <Isotope>s and <Material>s.
!=====================================================================
subroutine read_material(index, words, n_words)
integer, intent(in) :: index
character(*), intent(in) :: words(n_words)
integer, intent(in) :: n_words
character(100) :: msg
integer :: i
integer :: n_isotopes
type(Material), pointer :: mat
if (mod(n_words,2) /= 0) then
msg = "Invalid number of arguments for material: " // words(2)
call error(msg)
end if
n_isotopes = (n_words-2)/2
mat => materials(index)
mat%n_isotopes = n_isotopes
! allocate isotope and density list
allocate(mat%names(n_isotopes))
allocate(mat%isotopes(n_isotopes))
allocate(mat%table(n_isotopes))
allocate(mat%atom_percent(n_isotopes))
! 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))
end do
end subroutine read_material
!=====================================================================
! NORMALIZE_AO normalizes the atom or weight percentages for each
! material
!=====================================================================
subroutine normalize_ao()
type(xsData), pointer :: iso
type(Material), pointer :: mat
integer :: index
integer :: i, j
character(10) :: key
real(8) :: sum_percent
real(8) :: awr, w
character(100) :: msg
! first find the index in the xsdata array for each isotope in
! each material
do i = 1, n_materials
mat => materials(i)
! Check to make sure either all atom percents or all weight
! percents are given
if (.not. (all(mat%atom_percent > 0) .or. &
& all(mat%atom_percent < 0))) then
msg = "Cannot mix atom and weight percents in material " // &
& int_to_str(mat%uid)
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
! 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
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)
awr = xsdatas(index)%awr
w = mat%atom_percent(j)
mat%atom_percent(j) = (w/awr)/sum_percent
end do
end do
end subroutine normalize_ao
!=====================================================================
! READ_XS_LIBRARY parses the data on a xs_library card. This card
! specifies what cross section library should be used by default
@ -373,29 +508,121 @@ contains
! Read number of cycles
read(words(2), fmt='(I8)', iostat=readError) n_cycles
if ( readError > 0 ) then
if (readError > 0) then
msg = "Invalid number of cycles: " // words(2)
call error( msg )
call error(msg)
end if
! Read number of inactive cycles
read(words(3), fmt='(I8)', iostat=readError) n_inactive
if ( readError > 0 ) then
if (readError > 0) then
msg = "Invalid number of inactive cycles: " // words(2)
call error( msg )
call error(msg)
end if
! Read number of particles
read(words(4), fmt='(I8)', iostat=readError) n_particles
if ( readError > 0 ) then
msg = "Invalid number of particles: " // words(2)
call error( msg )
call error(msg)
end if
end subroutine read_criticality
!=====================================================================
! READ_LINE reads a line from a file open on a unit
!=====================================================================
subroutine read_line(unit, line, ioError)
integer, intent(in) :: unit
character(max_line), intent(out) :: line
integer, intent(out) :: ioError
read(UNIT=unit, FMT='(A250)', IOSTAT=ioError) line
end subroutine read_line
!=====================================================================
! SKIP_LINES skips 'n_lines' lines from a file open on a unit
!=====================================================================
subroutine skip_lines(unit, n_lines, ioError)
integer, intent(in) :: unit
integer, intent(in) :: n_lines
integer, intent(out) :: ioError
integer :: i
character(max_line) :: tmp
do i = 1, n_lines
read(UNIT=unit, fmt='(A250)', IOSTAT=ioError) tmp
end do
end subroutine skip_lines
!=====================================================================
! READ_DATA_INT reads integer data into an arrya from a file open
!=====================================================================
subroutine read_data_int(in, array, n, lines, words_per_line)
integer, intent(in) :: in
integer, intent(out) :: array(n)
integer, intent(in) :: n
integer, intent(in) :: lines
integer, intent(in) :: words_per_line
integer :: i, j
character(250) :: line
character(32) :: words(max_words)
integer :: ioError
integer :: n_words
integer :: index
integer :: n_last
index = 1
do i = 1, lines
if (i == lines) then
read(in,*) array(index:n)
else
read(in,*) array(index:index+words_per_line-1)
index = index + words_per_line
end if
end do
end subroutine read_data_int
!=====================================================================
! READ_DATA_REAL reads real(8) data into an arrya from a file open
!=====================================================================
subroutine read_data_real(in, array, n, lines, words_per_line)
integer, intent(in) :: in
real(8), intent(out) :: array(n)
integer, intent(in) :: n
integer, intent(in) :: lines
integer, intent(in) :: words_per_line
integer :: i, j
character(250) :: line
character(32) :: words(max_words)
integer :: ioError
integer :: n_words
integer :: index
index = 1
do i = 1, lines
if (i == lines) then
read(in,*) array(index:n)
else
read(in,*) array(index:index+words_per_line-1)
index = index + words_per_line
end if
end do
end subroutine read_data_real
end module fileio

97
src/global.f90
View file

@ -5,13 +5,28 @@ module global
implicit none
! Main arrays for cells, surfaces, materials
type(Cell), allocatable, target :: cells(:)
type(Surface), allocatable, target :: surfaces(:)
type(Material), allocatable, target :: materials(:)
type(Cell), allocatable, target :: cells(:)
type(Surface), allocatable, target :: surfaces(:)
type(Material), allocatable, target :: materials(:)
type(Isotope), allocatable, target :: isotopes(:)
type(xsData), allocatable, target :: xsdatas(:)
integer :: n_cells ! # of cells
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
! Cross section arrays
type(AceContinuous), pointer :: xs_continuous(:)
type(AceThermal), allocatable, target :: xs_thermal(:)
! unionized energy grid
real(8), allocatable :: energy_grid(:)
! Histories/cycles/etc for both external source and criticality
integer :: n_particles ! # of particles (per cycle for criticality)
integer :: n_cycles ! # of cycles
@ -24,10 +39,19 @@ module global
type(Neutron), allocatable, target :: source_bank(:)
type(Bank), allocatable, target :: fission_bank(:)
! Paths to input file, cross section data, etc
character(100) :: &
& path_input, &
& path_xsdata
! Physical constants
real(8), parameter :: &
& PI = 2.*acos(0.0), & ! pi
& INFINITY = huge(0.0_8) ! positive infinity
real(8), parameter :: &
& PI = 2.*acos(0.0), & ! pi
& MASS_NEUTRON = 1.0086649156, & ! mass of a neutron
& MASS_PROTON = 1.00727646677, & ! mass of a proton
& AMU = 1.66053873e-27, & ! 1 amu in kg
& K_BOLTZMANN = 8.617342e-5, & ! Boltzmann constant in eV/K
& INFINITY = huge(0.0_8) ! positive infinity
! Boundary conditions
integer, parameter :: &
@ -81,18 +105,14 @@ module global
& PHOTON_ = 2, &
& ELECTRON_ = 3
character(32) :: inputfile
integer :: verbosity = 5
integer, parameter :: max_words = 100
integer, parameter :: max_line = 250
! Versioning numbers
integer, parameter :: VERSION_MAJOR = 0
integer, parameter :: VERSION_MINOR = 1
integer, parameter :: VERSION_RELEASE = 1
! Key length for dictionary
integer, parameter :: DICT_KEY_LENGTH = 20
integer, parameter :: VERSION_RELEASE = 2
contains
@ -118,23 +138,16 @@ contains
subroutine free_memory()
! Deallocate cells, surfaces, materials
if (allocated(cells)) then
deallocate(cells)
end if
if (allocated(surfaces)) then
deallocate(surfaces)
end if
if (allocated(materials)) then
deallocate(materials)
end if
if (allocated(cells)) deallocate(cells)
if (allocated(surfaces)) deallocate(surfaces)
if (allocated(materials)) deallocate(materials)
! Deallocate xsdata list
if (allocated(xsdatas)) deallocate(xsdatas)
! Deallocate fission and source bank
if (allocated(fission_bank)) then
deallocate(fission_bank)
end if
if (allocated(source_bank)) then
deallocate(source_bank)
end if
if (allocated(fission_bank)) deallocate(fission_bank)
if (allocated(source_bank)) deallocate(source_bank)
! End program
stop
@ -146,14 +159,36 @@ contains
! to integers less than 10 billion.
!=====================================================================
function int_to_str( num )
function int_to_str(num) result(str)
integer, intent(in) :: num
character(10) :: int_to_str
character(10) :: str
write ( int_to_str, '(I10)' ) num
int_to_str = adjustl(int_to_str)
write ( str, '(I10)' ) num
str = adjustl(str)
end function int_to_str
!=====================================================================
! INT_TO_STR converts a string to an integer.
!=====================================================================
function str_to_int(str) result(num)
character(*), intent(in) :: str
integer :: num
character(5) :: fmt
integer :: w
w = len_trim(str)
! Create format specifier for reading string
write(fmt, '("(I",I2,")")') w
! read string into integer
read(str,fmt) num
end function str_to_int
end module global

3
src/input_sample
View file

@ -6,9 +6,10 @@ cell 2 1 2 -1
surface 1 sph 0 0 0 5
surface 2 sph 0 0.0 0.0 3
material 1 U235 1.00
material 1 1001.03c 2.0 8016.03c 1.0 92235.12c 3.0
source box -3 -3 -3 3 3 3
xs_library endfb7
xs_data /opt/serpent/xsdata/endfb7
criticality 1 1 10

19
src/main.f90
View file

@ -1,7 +1,8 @@
program main
use global
use fileio, only: read_input, read_command_line, read_count
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 mcnp_random, only: RN_init_problem, rang, RN_init_particle
@ -9,6 +10,9 @@ program main
use physics, only: transport
use data_structures, only: Dictionary, dict_create, dict_add_key, &
& dict_get_key
use cross_section, only: read_xsdata
use ace, only: read_xs
use energy_grid, only: unionized_grid
implicit none
@ -17,6 +21,7 @@ program main
! Print the OpenMC title and version/date/time information
call title()
verbosity = 10
! Initialize random number generator
call RN_init_problem( 3, 0_8, 0_8, 0_8, 0 )
@ -30,15 +35,21 @@ program main
! Read input file -- make a first pass through the file to count
! cells, surfaces, etc in order to allocate arrays, then do a second
! pass to actually read values
call read_count(inputfile)
call read_input(inputfile)
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
call echo_input()
! create source particles
call init_source()
! start problem
verbosity = 10
surfaces(1)%bc = BC_VACUUM
call run_problem()

4
src/output.f90
View file

@ -203,13 +203,13 @@ module output
if ( val(2) < 10 ) then
if ( val(3) < 10 ) then
today_date = date(6:6) // "/" // date(7:7) // "/" // date(1:4)
today_date = date(6:6) // "/" // date(8:8) // "/" // date(1:4)
else
today_date = date(6:6) // "/" // date(7:8) // "/" // date(1:4)
end if
else
if ( val(3) < 10 ) then
today_date = date(5:6) // "/" // date(7:7) // "/" // date(1:4)
today_date = date(5:6) // "/" // date(8:8) // "/" // date(1:4)
else
today_date = date(5:6) // "/" // date(7:8) // "/" // date(1:4)
end if

54
src/search.f90 Normal file
View file

@ -0,0 +1,54 @@
module search
contains
!=====================================================================
! BINARY_SEARCH performs a binary search of an array to find where a
! specific value lies in the array. This is used extensively for
! energy grid searching
!=====================================================================
function binary_search( array, n, val ) result( index )
real(8), intent(in) :: array(n)
integer, intent(in) :: n
real(8), intent(in) :: val
integer :: index
integer :: L
integer :: R
real(8) :: testval
L = 1
R = n
if (val < array(L) .or. val > array(R)) then
! error
end if
do while (R - L > 1)
! Check boundaries
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
return
end if
! Find values at midpoint
index = L + (R - L)/2
testval = array(index)
if (val > testval) then
L = index + 1
elseif (val < testval) then
R = index - 1
end if
end do
index = L
end function binary_search
end module search

8
src/string.f90
View file

@ -1,7 +1,7 @@
module string
use global, only: max_words
use output, only: error
use output, only: error, warning
implicit none
@ -27,6 +27,7 @@ contains
integer :: i ! current index
integer :: i_start ! starting index of word
integer :: i_end ! ending index of word
character(250) :: msg
i_start = 0
i_end = 0
@ -43,6 +44,11 @@ contains
if (i == len(string)) i_end = i
if (i_end > 0) then
n = n + 1
if ( i_end - i_start + 1 > len(words(n)) ) then
msg = "The word '" // string(i_start:i_end) // "' is longer than " &
& // "the space allocated for it."
call warning( msg )
end if
words(n) = string(i_start:i_end)
! reset indices
i_start = 0

191
src/types.f90
View file

@ -2,7 +2,12 @@ module types
implicit none
type :: Surface
!=====================================================================
! SURFACE type defines a first- or second-order surface that can be
! used to construct closed volumes (cells)
!=====================================================================
type Surface
integer :: uid
integer :: type
real(8), allocatable :: coeffs(:)
@ -11,14 +16,23 @@ module types
integer :: bc
end type Surface
type :: Cell
!=====================================================================
! CELL defines a closed volume by its bounding surfaces
!=====================================================================
type Cell
integer :: uid
integer :: n_items
integer, allocatable :: boundary_list(:)
integer :: material
end type Cell
type :: Neutron
!=====================================================================
! NEUTRON describes the state of a neutron being transported through
! the geometry
!=====================================================================
type Neutron
integer :: uid ! Unique ID
real(8) :: xyz(3) ! location
real(8) :: uvw(3) ! directional cosines
@ -28,27 +42,188 @@ module types
logical :: alive ! is particle alive?
end type Neutron
type :: Bank
!=====================================================================
! BANK is used for storing fission sites in criticality
! calculations. Since all the state information of a neutron is not
! needed, this type allows sites to be stored with less memory
!=====================================================================
type Bank
integer :: uid ! Unique ID
real(8) :: xyz(3) ! Location of bank particle
end type Bank
type :: Isotope
!=====================================================================
! ISOTOPE describes an isotope, e.g. U-235, within a material. Note
! that two separate variables must be used for the same isotope in two
! different materials since they will generally have different
! densities
!=====================================================================
type Isotope
integer :: uid ! unique identifier
integer :: zaid ! ZAID, e.g. 92235
character(3) :: xs ! cross section identifier, e.g. 70c
real(8) :: density ! density in atom/b-cm
end type Isotope
type :: Material
!=====================================================================
! MATERIAL describes a material by its constituent isotopes
!=====================================================================
type Material
integer :: uid
integer :: n_isotopes
type(Isotope), allocatable :: isotopes(:)
character(10), allocatable :: names(:)
integer, allocatable :: isotopes(:)
integer, allocatable :: table(:)
real(8), allocatable :: atom_percent(:)
integer :: sab_table
end type Material
type :: ExtSource
!=====================================================================
! EXTSOURCE describes an external source of neutrons for a
! fixed-source problem or for the starting source in a criticality
! problem
!=====================================================================
type ExtSource
integer :: type ! type of source, e.g. 'box' or 'cell'
real(8), allocatable :: values(:) ! values for particular source type
end type ExtSource
!=====================================================================
! ACEREACTION contains the cross-section and secondary energy and
! angle distributions for a single reaction in a continuous-energy
! ACE-format table
!=====================================================================
type AceReaction
integer :: MT
real(8) :: Q_value
real(8) :: TY
integer :: energy_index
real(8), allocatable :: sigma(:)
real(8), allocatable :: ang_cos(:,:)
real(8), allocatable :: ang_pdf(:,:)
real(8), allocatable :: ang_cdf(:,:)
end type AceReaction
!=====================================================================
! ACECONTINUOUS contains all the data for an ACE-format
! continuous-energy cross section. The ACE format (A Compact ENDF
! format) is used in MCNP and several other Monte Carlo codes.
!=====================================================================
type AceContinuous
character(20) :: name
real(8) :: awr
real(8) :: temp
real(8), allocatable :: energy(:)
real(8), allocatable :: sigma_t(:)
real(8), allocatable :: sigma_a(:)
real(8), allocatable :: sigma_el(:)
real(8), allocatable :: heating(:)
integer :: nu_p_type
real(8), allocatable :: nu_p_energy(:)
real(8), allocatable :: nu_p_value(:)
integer :: nu_t_type
real(8), allocatable :: nu_t_energy(:)
real(8), allocatable :: nu_t_value(:)
real(8), allocatable :: nu_d_energy(:)
real(8), allocatable :: nu_d_value(:)
real(8), allocatable :: nu_d_precursor_const(:,:)
real(8), allocatable :: nu_d_precursor_energy(:,:)
real(8), allocatable :: nu_d_precursor_prob(:,:)
type(AceReaction), pointer :: reactions(:)
end type AceContinuous
!=====================================================================
! ACETHERMAL contains S(a,b) data for thermal neutron scattering,
! typically off of light isotopes such as water, graphite, Be, etc
!=====================================================================
type AceThermal
character(20) :: name
real(8) :: awr
real(8) :: temperature
real(8), allocatable :: inelastic_e_in(:)
real(8), allocatable :: inelastic_sigma(:)
real(8), allocatable :: inelastic_e_out(:,:)
real(8), allocatable :: inelastic_mu_out(:,:)
real(8), allocatable :: elastic_e_in(:)
real(8), allocatable :: elastic_P(:)
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
!=====================================================================
type xsData
character(10) :: alias
character(10) :: id
integer :: type
integer :: zaid
integer :: isomeric
real(8) :: awr
real(8) :: temp
integer :: binary
character(100) :: path
end type xsData
end module types

49
src/unittest.f90 Normal file
View file

@ -0,0 +1,49 @@
program unittest
use global
use energy_grid, only: add_grid_points
call test_energy_grid()
contains
!=====================================================================
! TEST_ENERGY_GRID is a test of the energy grid unionization
! subroutine on a simplified set of data
!=====================================================================
subroutine test_energy_grid()
type(LinkedListGrid), pointer :: list => null()
type(LinkedListGrid), pointer :: current => null()
real(8), allocatable :: energy(:)
integer :: n, i
integer :: n_list
n = 10
allocate(energy(n))
do i = 1,n
energy(i) = 1.1248_8*i
end do
! create list
allocate(list)
current => list
current%energy = 4.23_8
allocate(current%next)
current => current%next
current%energy = 7.23_8
call add_grid_points(list, energy, n)
current => list
do while(associated(current))
print *,current%energy
current => current%next
end do
deallocate(energy)
end subroutine test_energy_grid
end program unittest