crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-29 06:35:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Double-indexing for energy grid, made Lists and Dictionaries more generic with interfaces. Sampling of nuclides.

Browse source
This commit is contained in:
Paul Romano 2011-02-10 18:56:37 +00:00
parent e64b4f53b9
commit 2610260f56
16 changed files with 1711 additions and 603 deletions
Download patch file Download diff file Expand all files Collapse all files

31
ChangeLog
View file

@ -1,3 +1,34 @@
2011-02-09 Paul Romano <romano7@mit.edu>
* 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 <romano7@mit.edu>
* ace.f90: Added initial routines for reading continuous-energy

11
src/Makefile
View file

@ -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:

89
src/ace.f90
View file

@ -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

66
src/cross_section.f90
View file

@ -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

985
src/data_structures.f90
View file

File diff suppressed because it is too large Load diff

86
src/energy_grid.f90
View file

@ -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

275
src/fileio.f90
View file

@ -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
! <Isotope>s and <Material>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

323
src/geometry.f90
View file

@ -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

29
src/global.f90
View file

@ -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

14
src/input_sample
View file

@ -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

41
src/main.f90
View file

@ -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()

161
src/physics.f90
View file

@ -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

16
src/search.f90
View file

@ -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

4
src/source.f90
View file

@ -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)

6
src/string.f90
View file

@ -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

177
src/types.f90
View file

@ -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