mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 21:55:41 -04:00
Added initial universe capability (not complete), and fixed angular distribution sampling.
This commit is contained in:
Paul Romano
parent
59f08f8a08
commit
03a2634fb4
12 changed files with 645 additions and 74 deletions
23
ChangeLog
23
ChangeLog
|
|
@ -1,3 +1,26 @@
|
|||
2011-02-27 Paul Romano <romano7@mit.edu>
|
||||
|
||||
* ace.f90: Fixed minor bugs associated with reading angular
|
||||
distributions.
|
||||
* data_structures.f90: Added dict_keys routines.
|
||||
* fileio.f90: Considerably changed read_count and read_cell to be
|
||||
able to handle universes.
|
||||
* geometry.f90: Added universe logic to find_cell and allowed it
|
||||
to be recursively called.
|
||||
* global.f90: Added universe array and dictionary, cell types,
|
||||
pointer to lowest universe, angular distribution types,
|
||||
interpolation types, as well as ERROR_CODE and adjusted str_to_int
|
||||
to return ERROR_CODE on error.
|
||||
* output.f90: Added print_cell, print_universe, and print_summary
|
||||
subroutines. Need a few more useful routines like this!
|
||||
* physics.f90: Added ability to sample cosine of scattering angle
|
||||
from a 32-equiprobable bin distribution or a tabular
|
||||
distribution. This is done by the sample_energy routine. The
|
||||
rotate_angle subroutine takes a particles direction and applies
|
||||
that scattering angle to it. Also added some logic for universes.
|
||||
* types.f90: Changed some attributes in the Neutron structure,
|
||||
added type to angular distribution within AceReaction.
|
||||
|
||||
2011-02-22 Paul Romano <romano7@mit.edu>
|
||||
|
||||
* ace.f90: Added read_angular_dist and read_energy_dist
|
||||
|
|
|
|||
|
|
@ -59,7 +59,7 @@ 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
|
||||
output.o: global.o types.o data_structures.o
|
||||
physics.o: types.o global.o mcnp_random.o geometry.o output.o \
|
||||
search.o endf.o
|
||||
search.o: output.o
|
||||
|
|
|
|||
15
src/ace.f90
15
src/ace.f90
|
|
@ -365,6 +365,7 @@ contains
|
|||
NE = XSS(JXS9 + LOCB - 1)
|
||||
rxn % adist_n_energy = NE
|
||||
allocate(rxn % adist_energy(NE))
|
||||
allocate(rxn % adist_type(NE))
|
||||
allocate(rxn % adist_location(NE))
|
||||
|
||||
! read incoming energy grid and location of tables
|
||||
|
|
@ -376,11 +377,16 @@ contains
|
|||
size = 0
|
||||
do j = 1, NE
|
||||
LC = rxn % adist_location(j)
|
||||
if (LC > 0) then
|
||||
if (LC == 0) then
|
||||
! isotropic
|
||||
rxn % adist_type(j) = ANGLE_ISOTROPIC
|
||||
elseif (LC > 0) then
|
||||
! 32 equiprobable bins
|
||||
size = size + 32
|
||||
rxn % adist_type(j) = ANGLE_32_EQUI
|
||||
size = size + 33
|
||||
elseif (LC < 0) then
|
||||
! tabular distribution
|
||||
rxn % adist_type(j) = ANGLE_TABULAR
|
||||
NP = XSS(JXS9 + abs(LC))
|
||||
size = size + 2 + 3*NP
|
||||
end if
|
||||
|
|
@ -396,7 +402,10 @@ contains
|
|||
XSS_index = JXS9 + abs(LC) - 1
|
||||
rxn % adist_data = get_real(size)
|
||||
|
||||
! print *, trim(table%name), rxn%MT, size
|
||||
! change location pointers since they are currently relative to
|
||||
! JXS(9)
|
||||
LC = abs(rxn % adist_location(1))
|
||||
rxn % adist_location = abs(rxn % adist_location) - LC + 1
|
||||
|
||||
end do
|
||||
|
||||
|
|
|
|||
|
|
@ -18,8 +18,8 @@ module data_structures
|
|||
|
||||
implicit none
|
||||
|
||||
integer, parameter, private :: hash_size = 4993
|
||||
integer, parameter, private :: multiplier = 31
|
||||
integer, parameter :: hash_size = 4993
|
||||
integer, parameter :: multiplier = 31
|
||||
integer, parameter :: DICT_NULL = 0
|
||||
|
||||
!=====================================================================
|
||||
|
|
@ -82,6 +82,9 @@ module data_structures
|
|||
interface dict_hashkey
|
||||
module procedure dict_ci_hashkey, dict_ii_hashkey
|
||||
end interface
|
||||
interface dict_keys
|
||||
module procedure dict_ci_keys, dict_ii_keys
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
|
|
@ -1029,4 +1032,68 @@ contains
|
|||
|
||||
end function dict_ii_hashkey
|
||||
|
||||
!=====================================================================
|
||||
! DICT_CI_KEYS
|
||||
!=====================================================================
|
||||
|
||||
function dict_ci_keys(dict) result(head)
|
||||
|
||||
type(DictionaryCI), pointer :: dict
|
||||
type(ListKeyValueCI), pointer :: head
|
||||
type(ListKeyValueCI), pointer :: current => null()
|
||||
type(ListKeyValueCI), pointer :: elem => null()
|
||||
|
||||
integer :: i, j
|
||||
|
||||
head => null()
|
||||
|
||||
do i = 1, size(dict%table)
|
||||
elem => dict%table(i)%list
|
||||
do while (associated(elem))
|
||||
if (.not. associated(head)) then
|
||||
allocate(head)
|
||||
current => head
|
||||
else
|
||||
allocate(current%next)
|
||||
current => current%next
|
||||
end if
|
||||
current%data = elem%data
|
||||
elem => elem%next
|
||||
end do
|
||||
end do
|
||||
|
||||
end function dict_ci_keys
|
||||
|
||||
!=====================================================================
|
||||
! DICT_II_KEYS
|
||||
!=====================================================================
|
||||
|
||||
function dict_ii_keys(dict) result(head)
|
||||
|
||||
type(DictionaryII), pointer :: dict
|
||||
type(ListKeyValueII), pointer :: head
|
||||
type(ListKeyValueII), pointer :: current => null()
|
||||
type(ListKeyValueII), pointer :: elem => null()
|
||||
|
||||
integer :: i, j
|
||||
|
||||
head => null()
|
||||
|
||||
do i = 1, size(dict%table)
|
||||
elem => dict%table(i)%list
|
||||
do while (associated(elem))
|
||||
if (.not. associated(head)) then
|
||||
allocate(head)
|
||||
current => head
|
||||
else
|
||||
allocate(current%next)
|
||||
current => current%next
|
||||
end if
|
||||
current%data = elem%data
|
||||
elem => elem%next
|
||||
end do
|
||||
end do
|
||||
|
||||
end function dict_ii_keys
|
||||
|
||||
end module data_structures
|
||||
|
|
|
|||
141
src/fileio.f90
141
src/fileio.f90
|
|
@ -1,15 +1,22 @@
|
|||
module fileio
|
||||
|
||||
use global
|
||||
use types, only: Cell, Surface, ExtSource
|
||||
use types, only: Cell, Surface, ExtSource, ListKeyValueCI
|
||||
use string, only: split_string_wl, lower_case, str_to_real, split_string, &
|
||||
& concatenate
|
||||
use output, only: message, warning, error
|
||||
use data_structures, only: dict_create, dict_add_key, dict_get_key, &
|
||||
& DICT_NULL
|
||||
& dict_has_key, DICT_NULL, dict_keys, list_size
|
||||
|
||||
implicit none
|
||||
|
||||
type(DictionaryII), pointer :: & ! this dictionary is used to count
|
||||
& ucount_dict => null() ! the number of universes as well
|
||||
! as how many cells each universe
|
||||
! contains
|
||||
|
||||
integer, allocatable :: index_cell_in_univ(:)
|
||||
|
||||
!=====================================================================
|
||||
! 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
|
||||
|
|
@ -60,8 +67,8 @@ contains
|
|||
|
||||
!=====================================================================
|
||||
! READ_COUNT makes a first pass through the input file to determine
|
||||
! how many cells, surfaces, materials, sources, etc there are in the
|
||||
! problem.
|
||||
! how many cells, universes, surfaces, materials, sources, etc there
|
||||
! are in the problem.
|
||||
!=====================================================================
|
||||
|
||||
subroutine read_count(filename)
|
||||
|
|
@ -73,14 +80,23 @@ contains
|
|||
integer :: in = 7
|
||||
integer :: ioError
|
||||
integer :: n
|
||||
integer :: count
|
||||
integer :: universe_num
|
||||
type(ListKeyValueII), pointer :: key_list
|
||||
type(Universe), pointer :: univ
|
||||
integer :: index
|
||||
|
||||
msg = "First pass through input file..."
|
||||
call message(msg, 5)
|
||||
|
||||
n_cells = 0
|
||||
n_universes = 0
|
||||
n_surfaces = 0
|
||||
n_materials = 0
|
||||
|
||||
call dict_create(ucount_dict)
|
||||
call dict_create(universe_dict)
|
||||
|
||||
open(FILE=filename, UNIT=in, STATUS='old', &
|
||||
& ACTION='read', IOSTAT=ioError)
|
||||
if (ioError /= 0) then
|
||||
|
|
@ -96,6 +112,20 @@ contains
|
|||
select case (trim(words(1)))
|
||||
case ('cell')
|
||||
n_cells = n_cells + 1
|
||||
|
||||
! For cells, we also need to check if there's a new universe
|
||||
! -- also for every cell add 1 to the count of cells for the
|
||||
! specified universe
|
||||
universe_num = str_to_int(words(3))
|
||||
if (.not. dict_has_key(ucount_dict, universe_num)) then
|
||||
n_universes = n_universes + 1
|
||||
count = 1
|
||||
call dict_add_key(universe_dict, universe_num, n_universes)
|
||||
else
|
||||
count = 1 + dict_get_key(ucount_dict, universe_num)
|
||||
end if
|
||||
call dict_add_key(ucount_dict, universe_num, count)
|
||||
|
||||
case ('surface')
|
||||
n_surfaces = n_surfaces + 1
|
||||
case ('material')
|
||||
|
|
@ -123,14 +153,47 @@ contains
|
|||
|
||||
! Allocate arrays for cells, surfaces, etc.
|
||||
allocate(cells(n_cells))
|
||||
allocate(universes(n_universes))
|
||||
allocate(surfaces(n_surfaces))
|
||||
allocate(materials(n_materials))
|
||||
|
||||
! Also allocate a list for keeping track of where cells have been
|
||||
! assigned in each universe
|
||||
allocate(index_cell_in_univ(n_universes))
|
||||
index_cell_in_univ = 0
|
||||
|
||||
! Set up dictionaries
|
||||
call dict_create(cell_dict)
|
||||
call dict_create(surface_dict)
|
||||
call dict_create(material_dict)
|
||||
|
||||
! We also need to allocate the cell count lists for each
|
||||
! universe. The logic for this is a little more convoluted. In
|
||||
! universe_dict, the (key,value) pairs are the uid of the universe
|
||||
! and the index in the array. In ucount_dict, it's the uid of the
|
||||
! universe and the number of cells.
|
||||
|
||||
key_list => dict_keys(universe_dict)
|
||||
do while (associated(key_list))
|
||||
! find index of universe in universes array
|
||||
index = key_list%data%value
|
||||
univ => universes(index)
|
||||
univ % uid = key_list%data%key
|
||||
|
||||
! check for lowest level universe
|
||||
if (univ % uid == 0) BASE_UNIVERSE = index
|
||||
|
||||
! find cell count for this universe
|
||||
count = dict_get_key(ucount_dict, key_list%data%key)
|
||||
|
||||
! allocate cell list for universe
|
||||
allocate(univ % cells(count))
|
||||
univ % n_cells = count
|
||||
|
||||
! move to next universe
|
||||
key_list => key_list%next
|
||||
end do
|
||||
|
||||
end subroutine read_count
|
||||
|
||||
!=====================================================================
|
||||
|
|
@ -149,6 +212,7 @@ contains
|
|||
integer :: ioError
|
||||
integer :: n, i
|
||||
integer :: index_cell
|
||||
integer :: index_universe
|
||||
integer :: index_surface
|
||||
integer :: index_material
|
||||
integer :: index_source
|
||||
|
|
@ -157,6 +221,7 @@ contains
|
|||
call message(msg, 5)
|
||||
|
||||
index_cell = 0
|
||||
index_universe = 0
|
||||
index_surface = 0
|
||||
index_material = 0
|
||||
index_source = 0
|
||||
|
|
@ -205,6 +270,9 @@ contains
|
|||
|
||||
end do
|
||||
|
||||
! deallocate array
|
||||
deallocate(index_cell_in_univ)
|
||||
|
||||
close(UNIT=in)
|
||||
|
||||
call adjust_indices()
|
||||
|
|
@ -245,13 +313,17 @@ contains
|
|||
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)
|
||||
if (c % fill /= 0) then
|
||||
! cell is filled with universe -- do nothing
|
||||
else
|
||||
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 if
|
||||
c%material = index
|
||||
end do
|
||||
|
||||
end subroutine adjust_indices
|
||||
|
|
@ -268,32 +340,63 @@ contains
|
|||
|
||||
integer :: ioError
|
||||
integer :: i
|
||||
integer :: universe_num
|
||||
integer :: n_items
|
||||
character(250) :: msg
|
||||
character(32) :: word
|
||||
type(Cell), pointer :: this_cell => null()
|
||||
type(Universe), pointer :: this_univ => null()
|
||||
|
||||
this_cell => cells(index)
|
||||
|
||||
! Read cell identifier
|
||||
read(words(2), FMT='(I8)', IOSTAT=ioError) this_cell%uid
|
||||
if (ioError > 0) then
|
||||
this_cell % uid = str_to_int(words(2))
|
||||
if (this_cell % uid == ERROR_CODE) then
|
||||
msg = "Invalid cell name: " // words(2)
|
||||
call error(msg)
|
||||
end if
|
||||
call dict_add_key(cell_dict, this_cell%uid, index)
|
||||
|
||||
! Read cell universe
|
||||
universe_num = str_to_int(words(3))
|
||||
if (universe_num == ERROR_CODE) then
|
||||
msg = "Invalid universe: " // words(3)
|
||||
call error(msg)
|
||||
end if
|
||||
this_cell % universe = dict_get_key(universe_dict, universe_num)
|
||||
|
||||
! Read cell material
|
||||
read(words(3), FMT='(I8)', IOSTAT=ioError) this_cell%material
|
||||
if (trim(words(4)) == 'fill') then
|
||||
this_cell % type = CELL_FILL
|
||||
this_cell % material = 0
|
||||
n_items = n_words - 5
|
||||
|
||||
! find universe
|
||||
universe_num = str_to_int(words(5))
|
||||
if (universe_num == ERROR_CODE) then
|
||||
msg = "Invalid universe fill: " // words(5)
|
||||
call error(msg)
|
||||
end if
|
||||
this_cell % fill = dict_get_key(universe_dict, universe_num)
|
||||
|
||||
else
|
||||
this_cell % type = CELL_NORMAL
|
||||
this_cell % material = str_to_int(words(4))
|
||||
this_cell % fill = 0
|
||||
if (this_cell % material == ERROR_CODE) then
|
||||
msg = "Invalid material number: " // words(4)
|
||||
call error(msg)
|
||||
end if
|
||||
n_items = n_words - 4
|
||||
end if
|
||||
|
||||
! 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_items))
|
||||
do i = 1, n_items
|
||||
word = words(i+3)
|
||||
word = words(i+n_words-n_items)
|
||||
if (word(1:1) == '(') then
|
||||
this_cell%boundary_list(i) = OP_LEFT_PAREN
|
||||
elseif (word(1:1) == ')') then
|
||||
|
|
@ -303,10 +406,16 @@ contains
|
|||
elseif (word(1:1) == '#') then
|
||||
this_cell%boundary_list(i) = OP_DIFFERENCE
|
||||
else
|
||||
read(word, FMT='(I8)', IOSTAT=ioError) this_cell%boundary_list(i)
|
||||
this_cell%boundary_list(i) = str_to_int(word)
|
||||
end if
|
||||
end do
|
||||
|
||||
! Add cell to the cell list in the corresponding universe
|
||||
i = this_cell % universe
|
||||
this_univ => universes(i)
|
||||
index_cell_in_univ(i) = index_cell_in_univ(i) + 1
|
||||
this_univ % cells(index_cell_in_univ(i)) = index
|
||||
|
||||
end subroutine read_cell
|
||||
|
||||
!=====================================================================
|
||||
|
|
|
|||
|
|
@ -85,35 +85,47 @@ contains
|
|||
! FIND_CELL determines what cell a source neutron is in
|
||||
!=====================================================================
|
||||
|
||||
subroutine find_cell(neut)
|
||||
recursive subroutine find_cell(univ, neut, found)
|
||||
|
||||
type(Neutron), pointer :: neut
|
||||
type(Universe), pointer :: univ
|
||||
type(Neutron), pointer :: neut
|
||||
logical, intent(inout) :: found
|
||||
|
||||
type(Cell), pointer :: this_cell
|
||||
logical :: found_cell
|
||||
character(250) :: msg
|
||||
integer :: i
|
||||
character(250) :: msg ! error message
|
||||
integer :: i ! index over cells
|
||||
type(Cell), pointer :: this_cell !
|
||||
type(Universe), pointer :: lower_univ
|
||||
|
||||
found = .false.
|
||||
|
||||
! determine what region in
|
||||
do i = 1, n_cells
|
||||
this_cell => cells(i)
|
||||
do i = 1, univ % n_cells
|
||||
this_cell => cells(univ % cells(i))
|
||||
|
||||
if (cell_contains(this_cell, neut)) then
|
||||
neut%cell = i
|
||||
found_cell = .true.
|
||||
neut%cell = dict_get_key(cell_dict, this_cell % uid)
|
||||
|
||||
! set current pointers
|
||||
cCell => this_cell
|
||||
cMaterial => materials(cCell%material)
|
||||
exit
|
||||
! If this cell contains a universe of lattice, search for
|
||||
! the particle in that universe/lattice
|
||||
if (this_cell % fill > 0) then
|
||||
lower_univ => universes(this_cell % fill)
|
||||
call find_cell(lower_univ, neut, found)
|
||||
if (found) then
|
||||
exit
|
||||
else
|
||||
msg = "Could not locate neutron in universe: "
|
||||
call error(msg)
|
||||
end if
|
||||
else
|
||||
! set current pointers
|
||||
found = .true.
|
||||
cCell => this_cell
|
||||
cMaterial => materials(cCell%material)
|
||||
exit
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
|
||||
! if neutron couldn't be located, print error
|
||||
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)
|
||||
end if
|
||||
end do
|
||||
|
||||
end subroutine find_cell
|
||||
|
||||
|
|
|
|||
|
|
@ -6,16 +6,19 @@ module global
|
|||
|
||||
! Main arrays for cells, surfaces, materials
|
||||
type(Cell), allocatable, target :: cells(:)
|
||||
type(Universe), allocatable, target :: universes(:)
|
||||
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_universes ! # of universes
|
||||
integer :: n_surfaces ! # of surfaces
|
||||
integer :: n_materials ! # of materials
|
||||
|
||||
! These dictionaries provide a fast lookup mechanism
|
||||
type(DictionaryII), pointer :: cell_dict
|
||||
type(DictionaryII), pointer :: universe_dict
|
||||
type(DictionaryII), pointer :: surface_dict
|
||||
type(DictionaryII), pointer :: material_dict
|
||||
type(DictionaryII), pointer :: isotope_dict
|
||||
|
|
@ -30,6 +33,7 @@ module global
|
|||
|
||||
! Current cell, surface, material
|
||||
type(Cell), pointer :: cCell
|
||||
type(Universe), pointer :: cUniverse
|
||||
type(Surface), pointer :: cSurface
|
||||
type(Material), pointer :: cMaterial
|
||||
|
||||
|
|
@ -79,6 +83,15 @@ module global
|
|||
& OP_UNION = huge(0) - 2, & ! Union operator
|
||||
& OP_DIFFERENCE = huge(0) - 3 ! Difference operator
|
||||
|
||||
! Cell types
|
||||
integer, parameter :: &
|
||||
& CELL_NORMAL = 1, & ! Cell with a specified material
|
||||
& CELL_FILL = 2, & ! Cell filled by a separate universe
|
||||
& CELL_LATTICE = 3 ! Cell filled with a lattice
|
||||
|
||||
! array index of universe 0
|
||||
integer :: BASE_UNIVERSE
|
||||
|
||||
! Surface types
|
||||
integer, parameter :: &
|
||||
& SURF_PX = 1, & ! Plane parallel to x-plane
|
||||
|
|
@ -112,12 +125,28 @@ module global
|
|||
& PROB_SOURCE = 1, & ! External source problem
|
||||
& PROB_CRITICALITY = 2 ! Criticality problem
|
||||
|
||||
! Angular distribution type
|
||||
integer, parameter :: &
|
||||
& ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution
|
||||
& ANGLE_32_EQUI = 2, & ! 32 equiprobable bins
|
||||
& ANGLE_TABULAR = 3 ! Tabular angular distribution
|
||||
|
||||
! Interpolation flag
|
||||
integer, parameter :: &
|
||||
& HISTOGRAM = 1, & ! Histogram interpolation on angle
|
||||
& LINEARLINEAR = 2 ! Linear-linear interpolation on angle
|
||||
|
||||
! Particle type
|
||||
integer, parameter :: &
|
||||
& NEUTRON_ = 1, &
|
||||
& PHOTON_ = 2, &
|
||||
& ELECTRON_ = 3
|
||||
|
||||
! Integer code for read error -- better hope this number is never
|
||||
! used in an input file!
|
||||
integer, parameter :: &
|
||||
& ERROR_CODE = -huge(0)
|
||||
|
||||
integer :: verbosity = 5
|
||||
integer, parameter :: max_words = 500
|
||||
integer, parameter :: max_line = 250
|
||||
|
|
@ -192,15 +221,18 @@ contains
|
|||
integer :: num
|
||||
|
||||
character(5) :: fmt
|
||||
integer :: w
|
||||
integer :: w
|
||||
integer :: ioError
|
||||
|
||||
! Determine width of string
|
||||
w = len_trim(str)
|
||||
|
||||
! Create format specifier for reading string
|
||||
write(fmt, '("(I",I2,")")') w
|
||||
write(UNIT=fmt, FMT='("(I",I2,")")') w
|
||||
|
||||
! read string into integer
|
||||
read(str,fmt) num
|
||||
read(UNIT=str, FMT=fmt, IOSTAT=ioError) num
|
||||
if (ioError > 0) num = ERROR_CODE
|
||||
|
||||
end function str_to_int
|
||||
|
||||
|
|
|
|||
|
|
@ -1,10 +1,13 @@
|
|||
# test input file
|
||||
|
||||
cell 100 40 -1
|
||||
cell 200 41 1 -2
|
||||
cell 1 0 fill 37 -2
|
||||
cell 100 37 40 -1
|
||||
cell 101 37 41 1
|
||||
cell 2 0 40 2 -3
|
||||
|
||||
surface 1 sph 0 0 0 3
|
||||
surface 2 sph 0 0 0 5
|
||||
surface 3 sph 0 0 0 8
|
||||
|
||||
material 40 -4.5 &
|
||||
92238.03c 1.0
|
||||
|
|
@ -13,7 +16,7 @@ material 41 -1.0 &
|
|||
1001.03c 2.0 &
|
||||
8016.03c 1.0
|
||||
|
||||
source box -1 -1 -1 1 1 1
|
||||
source box -4 -4 -4 4 4 4
|
||||
|
||||
xs_library endfb7
|
||||
xs_data /opt/serpent/xsdata/endfb7
|
||||
|
|
|
|||
21
src/main.f90
21
src/main.f90
|
|
@ -1,17 +1,17 @@
|
|||
program main
|
||||
|
||||
use global
|
||||
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, 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: dict_create, dict_add_key, dict_get_key
|
||||
use fileio, only: read_input, read_command_line, read_count, &
|
||||
& normalize_ao
|
||||
use output, only: title, echo_input, message, warning, error, &
|
||||
& print_summary
|
||||
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 cross_section, only: read_xsdata, material_total_xs
|
||||
use ace, only: read_xs
|
||||
use energy_grid, only: unionized_grid, original_indices
|
||||
use ace, only: read_xs
|
||||
use energy_grid, only: unionized_grid, original_indices
|
||||
|
||||
implicit none
|
||||
|
||||
|
|
@ -60,6 +60,7 @@ program main
|
|||
call material_total_xs()
|
||||
|
||||
call echo_input()
|
||||
call print_summary()
|
||||
|
||||
! create source particles
|
||||
call init_source()
|
||||
|
|
|
|||
130
src/output.f90
130
src/output.f90
|
|
@ -2,6 +2,8 @@ module output
|
|||
|
||||
use ISO_FORTRAN_ENV
|
||||
use global
|
||||
use types, only: Cell, Universe, Surface
|
||||
use data_structures, only: dict_get_key
|
||||
|
||||
implicit none
|
||||
|
||||
|
|
@ -218,4 +220,132 @@ module output
|
|||
|
||||
end subroutine get_today
|
||||
|
||||
!=====================================================================
|
||||
! PRINT_CELL displays the attributes of a cell
|
||||
!=====================================================================
|
||||
|
||||
subroutine print_cell(c)
|
||||
|
||||
type(Cell), pointer :: c
|
||||
|
||||
integer :: ou
|
||||
integer :: temp
|
||||
integer :: i
|
||||
type(Universe), pointer :: u => null()
|
||||
type(Material), pointer :: m => null()
|
||||
character(250) :: string
|
||||
|
||||
ou = OUTPUT_UNIT
|
||||
|
||||
write(ou,*) 'Cell ' // int_to_str(c % uid)
|
||||
temp = dict_get_key(cell_dict, c % uid)
|
||||
write(ou,*) ' Array Index = ' // int_to_str(temp)
|
||||
u => universes(c % universe)
|
||||
write(ou,*) ' Universe = ' // int_to_str(u % uid)
|
||||
if (c % fill == 0) then
|
||||
write(ou,*) ' Fill = NONE'
|
||||
else
|
||||
u => universes(c % fill)
|
||||
write(ou,*) ' Fill = ' // int_to_str(u % uid)
|
||||
end if
|
||||
if (c % material == 0) then
|
||||
write(ou,*) ' Material = NONE'
|
||||
else
|
||||
m => materials(c % material)
|
||||
write(ou,*) ' Material = ' // int_to_str(m % uid)
|
||||
end if
|
||||
string = ""
|
||||
do i = 1, c % n_items
|
||||
select case (c % boundary_List(i))
|
||||
case (OP_LEFT_PAREN)
|
||||
string = trim(string) // ' ('
|
||||
case (OP_RIGHT_PAREN)
|
||||
string = trim(string) // ' )'
|
||||
case (OP_UNION)
|
||||
string = trim(string) // ' :'
|
||||
case (OP_DIFFERENCE)
|
||||
string = trim(string) // ' !'
|
||||
case default
|
||||
string = trim(string) // ' ' // int_to_str(c % boundary_list(i))
|
||||
end select
|
||||
end do
|
||||
write(ou,*) ' Surface Specification:' // trim(string)
|
||||
write(ou,*)
|
||||
|
||||
! nullify associated pointers
|
||||
nullify(u)
|
||||
nullify(m)
|
||||
|
||||
end subroutine print_cell
|
||||
|
||||
!=====================================================================
|
||||
! PRINT_UNIVERSE displays the attributes of a universe
|
||||
!=====================================================================
|
||||
|
||||
subroutine print_universe(univ)
|
||||
|
||||
type(Universe), pointer :: univ
|
||||
|
||||
integer :: ou
|
||||
integer :: i
|
||||
character(250) :: string
|
||||
type(Cell), pointer :: c
|
||||
|
||||
ou = OUTPUT_UNIT
|
||||
|
||||
write(ou,*) 'Universe ' // int_to_str(univ % uid)
|
||||
write(ou,*) ' Level = ' // int_to_str(univ % level)
|
||||
string = ""
|
||||
do i = 1, univ % n_cells
|
||||
c => cells(univ % cells(i))
|
||||
string = trim(string) // ' ' // int_to_str(c % uid)
|
||||
end do
|
||||
write(ou,*) ' Cells =' // trim(string)
|
||||
write(ou,*)
|
||||
|
||||
end subroutine print_universe
|
||||
|
||||
!=====================================================================
|
||||
! PRINT_SUMMARY displays the attributes of all cells, universes,
|
||||
! surfaces and materials read in the input file. Very useful for
|
||||
! debugging!
|
||||
!=====================================================================
|
||||
|
||||
subroutine print_summary()
|
||||
|
||||
type(Surface), pointer :: s
|
||||
type(Cell), pointer :: c
|
||||
type(Universe), pointer :: u
|
||||
type(Material), pointer :: m
|
||||
integer :: i
|
||||
integer :: ou
|
||||
|
||||
ou = OUTPUT_UNIT
|
||||
|
||||
! print summary of cells
|
||||
write(ou,*) '============================================='
|
||||
write(ou,*) '=> CELL SUMMARY <='
|
||||
write(ou,*) '============================================='
|
||||
write(ou,*)
|
||||
do i = 1, n_cells
|
||||
c => cells(i)
|
||||
call print_cell(c)
|
||||
end do
|
||||
|
||||
! print summary of universes
|
||||
write(ou,*) '============================================='
|
||||
write(ou,*) '=> UNIVERSE SUMMARY <='
|
||||
write(ou,*) '============================================='
|
||||
write(ou,*)
|
||||
do i = 1, n_universes
|
||||
u => universes(i)
|
||||
call print_universe(u)
|
||||
end do
|
||||
|
||||
! print summary of surfaces
|
||||
|
||||
! print summary of materials
|
||||
|
||||
end subroutine print_summary
|
||||
|
||||
end module output
|
||||
|
|
|
|||
173
src/physics.f90
173
src/physics.f90
|
|
@ -4,7 +4,7 @@ module physics
|
|||
use geometry, only: find_cell, dist_to_boundary, cross_boundary
|
||||
use types, only: Neutron
|
||||
use mcnp_random, only: rang
|
||||
use output, only: error, message
|
||||
use output, only: error, message, print_universe
|
||||
use search, only: binary_search
|
||||
use endf, only: reaction_name
|
||||
|
||||
|
|
@ -19,11 +19,12 @@ contains
|
|||
|
||||
subroutine transport(neut)
|
||||
|
||||
type(Neutron), pointer, intent(inout) :: neut
|
||||
type(Neutron), pointer :: neut
|
||||
|
||||
character(250) :: msg
|
||||
integer :: i, surf
|
||||
logical :: alive
|
||||
logical :: found_cell
|
||||
|
||||
real(8) :: d_to_boundary
|
||||
real(8) :: d_to_collision
|
||||
|
|
@ -32,11 +33,23 @@ contains
|
|||
real(8) :: f ! interpolation factor
|
||||
real(8) :: tmp(3)
|
||||
integer :: IE ! index on energy grid
|
||||
type(Universe), pointer :: univ
|
||||
|
||||
! determine what cell the particle is in
|
||||
if (neut%cell == 0) then
|
||||
call find_cell(neut)
|
||||
univ => universes(BASE_UNIVERSE)
|
||||
call find_cell(univ, neut, found_cell)
|
||||
print *, sqrt(neut%xyz(1)**2 + neut%xyz(2)**2 + neut%xyz(3)**2), cells(neut%cell) % uid
|
||||
|
||||
! if neutron couldn't be located, print error
|
||||
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)
|
||||
end if
|
||||
end if
|
||||
|
||||
return
|
||||
|
||||
if (verbosity >= 10) then
|
||||
msg = "=== Particle " // trim(int_to_str(neut%uid)) // " ==="
|
||||
call message(msg, 10)
|
||||
|
|
@ -268,6 +281,158 @@ contains
|
|||
|
||||
end subroutine n_gamma
|
||||
|
||||
!=====================================================================
|
||||
! SAMPLE_ANGLE samples the cosine of the angle between incident and
|
||||
! exiting particle directions either from 32 equiprobable bins or from
|
||||
! a tabular distribution.
|
||||
!=====================================================================
|
||||
|
||||
subroutine sample_angle(rxn, E)
|
||||
|
||||
type(AceReaction), pointer :: rxn ! reaction
|
||||
real(8), intent(in) :: E ! incoming energy
|
||||
|
||||
real(8) :: xi ! random number on [0,1)
|
||||
integer :: interp ! type of interpolation
|
||||
integer :: type ! angular distribution type
|
||||
integer :: i ! incoming energy bin
|
||||
integer :: n ! number of incoming energy bins
|
||||
integer :: loc ! location in data array
|
||||
integer :: np ! number of points in cos distribution
|
||||
integer :: bin ! cosine bin
|
||||
real(8) :: f ! interpolation factor
|
||||
real(8) :: mu0 ! cosine in bin b
|
||||
real(8) :: mu1 ! cosine in bin b+1
|
||||
real(8) :: mu ! final cosine sampled
|
||||
real(8) :: c ! cumulative distribution frequency
|
||||
real(8) :: p0,p1 ! probability distribution
|
||||
character(250) :: msg ! error message
|
||||
|
||||
! determine number of incoming energies
|
||||
n = rxn % adist_n_energy
|
||||
|
||||
! find energy bin and calculate interpolation factor -- if the
|
||||
! energy is outside the range of the tabulated energies, choose
|
||||
! the first or last bins
|
||||
if (E < rxn%adist_energy(1)) then
|
||||
i = 1
|
||||
f = 0.0
|
||||
elseif (E > rxn%adist_energy(n)) then
|
||||
i = n - 1
|
||||
f = 1.0
|
||||
else
|
||||
i = binary_search(rxn%adist_energy, n, E)
|
||||
f = (E - rxn % adist_energy(i)) / &
|
||||
& (rxn % adist_energy(i+1) - rxn % adist_energy(i))
|
||||
end if
|
||||
|
||||
! Sample between the ith and (i+1)th bin
|
||||
if (f > rang()) i = i + 1
|
||||
|
||||
! check whether this is a 32-equiprobable bin or a tabular
|
||||
! distribution
|
||||
loc = rxn % adist_location(i)
|
||||
type = rxn % adist_type(i)
|
||||
if (type == ANGLE_ISOTROPIC) then
|
||||
mu = 2.0_8 * rang() - 1
|
||||
elseif (type == ANGLE_32_EQUI) then
|
||||
! sample cosine bin
|
||||
xi = rang()
|
||||
bin = 1 + int(32.0_8*xi)
|
||||
|
||||
! calculate cosine
|
||||
mu0 = rxn % adist_data(loc + bin - 1)
|
||||
mu1 = rxn % adist_data(loc + bin)
|
||||
mu = mu0 + (32.0_8 * xi - bin) * (mu1 - mu0)
|
||||
|
||||
elseif (type == ANGLE_TABULAR) then
|
||||
interp = rxn % adist_data(loc)
|
||||
np = rxn % adist_data(loc+1)
|
||||
|
||||
! determine outgoing cosine bin
|
||||
xi = rang()
|
||||
do bin = loc+2, loc+1+NP
|
||||
c = rxn % adist_data(bin+2*np)
|
||||
if (xi > c) exit
|
||||
end do
|
||||
|
||||
p0 = rxn % adist_data(bin+np)
|
||||
mu0 = rxn % adist_data(bin)
|
||||
if (interp == HISTOGRAM) then
|
||||
! Histogram interpolation
|
||||
mu = mu0 + (xi - c)/p0
|
||||
|
||||
elseif (interp == LINEARLINEAR) then
|
||||
! Linear-linear interpolation -- not sure how you come about
|
||||
! the formula given in the MCNP manual
|
||||
p1 = rxn % adist_data(bin+np+1)
|
||||
mu1 = rxn % adist_data(bin+1)
|
||||
|
||||
f = (p1 - p0)/(mu1 - mu0)
|
||||
if (f == ZERO) then
|
||||
mu = mu0 + (xi - c)/p0
|
||||
else
|
||||
mu = mu0 + (sqrt(p0*p0 + 2*f*(xi - c))-p0)/f
|
||||
end if
|
||||
else
|
||||
msg = "Unknown interpolation type: " // trim(int_to_str(interp))
|
||||
call error(msg)
|
||||
end if
|
||||
|
||||
else
|
||||
msg = "Unknown angular distribution type: " // trim(int_to_str(type))
|
||||
call error(msg)
|
||||
end if
|
||||
|
||||
end subroutine sample_angle
|
||||
|
||||
!=====================================================================
|
||||
! ROTATE_ANGLE rotates direction cosines through a polar angle whose
|
||||
! cosine is mu and through an azimuthal angle sampled uniformly. Note
|
||||
! that this is done with direct sampling rather than rejection as is
|
||||
! done in MCNP and SERPENT.
|
||||
!=====================================================================
|
||||
|
||||
subroutine rotate_angle(mu, u, v, w)
|
||||
|
||||
real(8), intent(in) :: mu ! cosine of angle
|
||||
real(8), intent(inout) :: u
|
||||
real(8), intent(inout) :: v
|
||||
real(8), intent(inout) :: w
|
||||
|
||||
real(8) :: phi, sinphi, cosphi
|
||||
real(8) :: a,b
|
||||
real(8) :: u0, v0, w0
|
||||
|
||||
! Copy original directional cosines
|
||||
u0 = u
|
||||
v0 = v
|
||||
w0 = w
|
||||
|
||||
! Sample azimuthal angle in [0,2pi)
|
||||
phi = 2.0_8 * PI * rang()
|
||||
|
||||
! Precompute factors to save flops
|
||||
sinphi = sin(phi)
|
||||
cosphi = cos(phi)
|
||||
a = sqrt(1.0_8 - mu*mu)
|
||||
b = sqrt(1.0_8 - w*w)
|
||||
|
||||
! Need to treat special case where sqrt(1 - w**2) is close to zero
|
||||
! by expanding about the v component rather than the w component
|
||||
if (b > 1e-10) then
|
||||
u = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b
|
||||
v = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b
|
||||
w = mu*w0 - a*b*cosphi
|
||||
else
|
||||
b = sqrt(1.0_8 - v*v)
|
||||
u = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b
|
||||
v = mu*v0 - a*b*cosphi
|
||||
w = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b
|
||||
end if
|
||||
|
||||
end subroutine rotate_angle
|
||||
|
||||
!=====================================================================
|
||||
! MAXWELL_SPECTRUM samples an energy from the Maxwell fission
|
||||
! distribution based on a rejection sampling scheme. This is described
|
||||
|
|
|
|||
|
|
@ -2,6 +2,21 @@ module types
|
|||
|
||||
implicit none
|
||||
|
||||
!=====================================================================
|
||||
! UNIVERSE
|
||||
!=====================================================================
|
||||
|
||||
type Universe
|
||||
integer :: uid
|
||||
integer :: type
|
||||
integer :: level
|
||||
integer :: n_cells
|
||||
integer, allocatable :: cells(:)
|
||||
real(8) :: x0
|
||||
real(8) :: y0
|
||||
real(8) :: z0
|
||||
end type Universe
|
||||
|
||||
!=====================================================================
|
||||
! SURFACE type defines a first- or second-order surface that can be
|
||||
! used to construct closed volumes (cells)
|
||||
|
|
@ -22,9 +37,12 @@ module types
|
|||
|
||||
type Cell
|
||||
integer :: uid
|
||||
integer :: type
|
||||
integer :: universe
|
||||
integer :: fill
|
||||
integer :: material
|
||||
integer :: n_items
|
||||
integer, allocatable :: boundary_list(:)
|
||||
integer :: material
|
||||
end type Cell
|
||||
|
||||
!=====================================================================
|
||||
|
|
@ -33,16 +51,17 @@ module types
|
|||
!=====================================================================
|
||||
|
||||
type Neutron
|
||||
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?
|
||||
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 :: universe ! current universe
|
||||
integer :: surface ! current surface
|
||||
real(8) :: wgt ! particle weight
|
||||
logical :: alive ! is particle alive?
|
||||
end type Neutron
|
||||
|
||||
!=====================================================================
|
||||
|
|
@ -115,6 +134,7 @@ module types
|
|||
! Secondary angle distribution
|
||||
integer :: adist_n_energy ! # of incoming energies
|
||||
real(8), allocatable :: adist_energy(:) ! incoming energy grid
|
||||
integer, allocatable :: adist_type(:) ! type of distribution
|
||||
integer, allocatable :: adist_location(:) ! location of each table
|
||||
real(8), allocatable :: adist_data(:) ! angular distribution data
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue