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

Changed attribute names in Cell, Neutron is now Particle.

Browse source
This commit is contained in:
Paul Romano 2011-03-10 22:37:17 +00:00
parent 4d8cca60ec
commit 1c42639ef5
9 changed files with 286 additions and 262 deletions
Download patch file Download diff file Expand all files Collapse all files

17
ChangeLog
View file

@ -1,3 +1,20 @@
2011-03-10 Paul Romano <romano7@mit.edu>
* types.f90: The major changes are that the n_items and
boundary_list attributes of 'Cell' are now n_surfaces and
surfaces. Even though this doesn't quite make sense since
parentheses and union operators are not surfaces per-se, it is
easier to understand within other modules with these names. The
'Neutron' type is now called 'Particle' and will be used for all
different particles (eventually).
* fileio.f90: Changed aforementioned names.
* geometry.f90: Changed aforementioned names.
* global.f90: Changed aforementioned names.
* main.f90: Changed aforementioned names.
* output.f90: Changed aforementioned names.
* physics.f90: Changed aforementioned names.
* source.f90: Changed aforementioned names.
2011-03-04 Paul Romano <romano7@mit.edu>
* fileio.f90: Added build_universe routine which sets pointers

64
src/fileio.f90
View file

@ -305,8 +305,8 @@ contains
do i = 1, n_cells
! adjust boundary list
c => cells(i)
do j = 1, c%n_items
surf_num = c%boundary_list(j)
do j = 1, c%n_surfaces
surf_num = c%surfaces(j)
if (surf_num < OP_DIFFERENCE) then
index = dict_get_key(surface_dict, abs(surf_num))
if (index == DICT_NULL) then
@ -315,7 +315,7 @@ contains
& " specified on cell " // trim(int_to_str(c%uid))
call error(msg)
end if
c%boundary_list(j) = sign(index, surf_num)
c%surfaces(j) = sign(index, surf_num)
end if
end do
@ -383,21 +383,20 @@ contains
integer :: ioError
integer :: i
integer :: universe_num
integer :: n_items
integer :: n_surfaces
character(250) :: msg
character(32) :: word
type(Cell), pointer :: this_cell => null()
type(Universe), pointer :: this_univ => null()
character(32) :: word
type(Cell), pointer :: c => null()
this_cell => cells(index)
c => cells(index)
! Read cell identifier
this_cell % uid = str_to_int(words(2))
if (this_cell % uid == ERROR_CODE) then
c % uid = str_to_int(words(2))
if (c % 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)
call dict_add_key(cell_dict, c%uid, index)
! Read cell universe
universe_num = str_to_int(words(3))
@ -405,13 +404,13 @@ contains
msg = "Invalid universe: " // words(3)
call error(msg)
end if
this_cell % universe = dict_get_key(universe_dict, universe_num)
c % universe = dict_get_key(universe_dict, universe_num)
! Read cell material
if (trim(words(4)) == 'fill') then
this_cell % type = CELL_FILL
this_cell % material = 0
n_items = n_words - 5
c % type = CELL_FILL
c % material = 0
n_surfaces = n_words - 5
! find universe
universe_num = str_to_int(words(5))
@ -419,44 +418,43 @@ contains
msg = "Invalid universe fill: " // words(5)
call error(msg)
end if
this_cell % fill = dict_get_key(universe_dict, universe_num)
c % 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
c % type = CELL_NORMAL
c % material = str_to_int(words(4))
c % fill = 0
if (c % material == ERROR_CODE) then
msg = "Invalid material number: " // words(4)
call error(msg)
end if
n_items = n_words - 4
n_surfaces = n_words - 4
end if
! Assign number of items
this_cell%n_items = n_items
c%n_surfaces = n_surfaces
! Read list of surfaces
allocate(this_cell%boundary_list(n_items))
do i = 1, n_items
word = words(i+n_words-n_items)
allocate(c%surfaces(n_surfaces))
do i = 1, n_surfaces
word = words(i+n_words-n_surfaces)
if (word(1:1) == '(') then
this_cell%boundary_list(i) = OP_LEFT_PAREN
c % surfaces(i) = OP_LEFT_PAREN
elseif (word(1:1) == ')') then
this_cell%boundary_list(i) = OP_RIGHT_PAREN
c % surfaces(i) = OP_RIGHT_PAREN
elseif (word(1:1) == ':') then
this_cell%boundary_list(i) = OP_UNION
c % surfaces(i) = OP_UNION
elseif (word(1:1) == '#') then
this_cell%boundary_list(i) = OP_DIFFERENCE
c % surfaces(i) = OP_DIFFERENCE
else
this_cell%boundary_list(i) = str_to_int(word)
c % surfaces(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)
i = c % universe
index_cell_in_univ(i) = index_cell_in_univ(i) + 1
this_univ % cells(index_cell_in_univ(i)) = index
universes(i) % cells(index_cell_in_univ(i)) = index
end subroutine read_cell

260
src/geometry.f90
View file

@ -13,28 +13,28 @@ contains
! CELL_CONTAINS determines whether a given point is inside a cell
!=====================================================================
function cell_contains(c, neut) result(in_cell)
function cell_contains(c, p) result(in_cell)
type(Cell), pointer :: c
type(Neutron), pointer :: neut
logical :: in_cell
type(Cell), pointer :: c
type(Particle), pointer :: p
logical :: in_cell
integer, allocatable :: expression(:)
integer :: specified_sense
integer :: actual_sense
integer :: n_boundaries
integer :: n_surfaces
integer :: i, j
integer :: surf_num
integer :: current_surface
type(Surface), pointer :: surf => null()
character(250) :: msg
current_surface = neut%surface
current_surface = p%surface
n_boundaries = size(c%boundary_list)
allocate(expression(n_boundaries))
expression = c%boundary_list
do i = 1, n_boundaries
n_surfaces = size(c%surfaces)
allocate(expression(n_surfaces))
expression = c%surfaces
do i = 1, n_surfaces
! Don't change logical operator
if (expression(i) >= OP_DIFFERENCE) then
@ -47,11 +47,11 @@ contains
! Check if the particle is currently on the specified surface
if (surf_num == current_surface) then
! neutron is on specified surface heading into cell
! particle 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
! particle is on specified surface, but heading other
! direction
expression(i) = 0
cycle
@ -59,7 +59,7 @@ contains
! Compare sense of point to specified sense
specified_sense = sign(1,expression(i))
call sense(surf, neut%xyz, actual_sense)
call sense(surf, p%xyz, actual_sense)
if (actual_sense == specified_sense) then
expression(i) = 1
else
@ -82,44 +82,48 @@ contains
end function cell_contains
!=====================================================================
! FIND_CELL determines what cell a source neutron is in
! FIND_CELL determines what cell a source particle is in within a
! particular universe. If the base universe is passed, the particle
! should be found as long as it's within the geometry
!=====================================================================
recursive subroutine find_cell(univ, neut, found)
recursive subroutine find_cell(univ, p, found)
type(Universe), pointer :: univ
type(Neutron), pointer :: neut
logical, intent(inout) :: found
type(Universe), pointer :: univ ! universe to search in
type(Particle), pointer :: p ! pointer to particle
logical, intent(inout) :: found ! particle found?
character(250) :: msg ! error message
integer :: i ! index over cells
type(Cell), pointer :: this_cell !
type(Universe), pointer :: lower_univ
type(Cell), pointer :: c ! pointer to cell
type(Universe), pointer :: lower_univ ! if particle is in lower
! universe, use this pointer
! to call recursively
found = .false.
! determine what region in
do i = 1, univ % n_cells
this_cell => cells(univ % cells(i))
c => cells(univ % cells(i))
if (cell_contains(this_cell, neut)) then
neut%cell = dict_get_key(cell_dict, this_cell % uid)
if (cell_contains(c, p)) then
p%cell = dict_get_key(cell_dict, c % uid)
! 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 (c % fill > 0) then
lower_univ => universes(c % fill)
call find_cell(lower_univ, p, found)
if (found) then
exit
else
msg = "Could not locate neutron in universe: "
msg = "Could not locate particle in universe: "
call error(msg)
end if
else
! set current pointers
found = .true.
cCell => this_cell
cCell => c
cMaterial => materials(cCell%material)
exit
end if
@ -130,12 +134,12 @@ contains
end subroutine find_cell
!=====================================================================
! CROSS_BOUNDARY moves a neutron into a new cell
! CROSS_BOUNDARY moves a particle into a new cell
!=====================================================================
subroutine cross_boundary(neut)
subroutine cross_boundary(p)
type(Neutron), pointer :: neut
type(Particle), pointer :: p
type(Surface), pointer :: surf
type(Cell), pointer :: c
@ -145,7 +149,7 @@ contains
logical :: found
type(Universe), pointer :: lower_univ => null()
surf => surfaces(abs(neut%surface))
surf => surfaces(abs(p%surface))
if (verbosity >= 10) then
msg = " Crossing surface " // trim(int_to_str(surf%uid))
call message(msg, 10)
@ -153,7 +157,7 @@ contains
! check for leakage
if (surf%bc == BC_VACUUM) then
neut%alive = .false.
p%alive = .false.
if (verbosity >= 10) then
msg = " Leaked out of surface " // trim(int_to_str(surf%uid))
call message(msg, 10)
@ -161,46 +165,46 @@ contains
return
end if
if (neut%surface > 0 .and. allocated(surf%neighbor_pos)) then
if (p%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)
index_cell = surf%neighbor_pos(i)
c => cells(index_cell)
if (cell_contains(c, neut)) then
if (cell_contains(c, p)) then
if (c % fill > 0) then
lower_univ => universes(c % fill)
call find_cell(lower_univ, neut, found)
call find_cell(lower_univ, p, found)
if (.not. found) then
msg = "Could not locate neutron in universe: "
msg = "Could not locate particle in universe: "
call error(msg)
end if
else
! set current pointers
neut%cell = index_cell
p%cell = index_cell
cCell => c
cMaterial => materials(cCell%material)
end if
return
end if
end do
elseif (neut%surface < 0 .and. allocated(surf%neighbor_neg)) then
elseif (p%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)
index_cell = surf%neighbor_neg(i)
c => cells(index_cell)
if (cell_contains(c, neut)) then
if (cell_contains(c, p)) then
if (c % fill > 0) then
lower_univ => universes(c % fill)
call find_cell(lower_univ, neut, found)
call find_cell(lower_univ, p, found)
if (.not. found) then
msg = "Could not locate neutron in universe: "
msg = "Could not locate particle in universe: "
call error(msg)
end if
else
! set current pointers
neut%cell = index_cell
p%cell = index_cell
cCell => c
cMaterial => materials(cCell%material)
end if
@ -213,8 +217,8 @@ contains
! cells
do i = 1, size(cells)
c => cells(i)
if (cell_contains(c, neut)) then
neut%cell = i
if (cell_contains(c, p)) then
p%cell = i
cCell => c
cMaterial => materials(cCell%material)
return
@ -222,7 +226,7 @@ contains
end do
! Couldn't find next cell anywhere!
msg = "After neutron crossed surface " // trim(int_to_str(neut%surface)) // &
msg = "After particle crossed surface " // trim(int_to_str(p%surface)) // &
& ", it could not be located in any cell and it did not leak."
call error(msg)
@ -230,14 +234,14 @@ contains
!=====================================================================
! DIST_TO_BOUNDARY calculates the distance to the nearest boundary of
! the cell 'cl' for a particle 'neut' traveling in a certain
! the cell 'cl' for a particle 'p' traveling in a certain
! direction. For a cell in a subuniverse that has a parent cell, also
! include the surfaces of the edge of the universe.
!=====================================================================
subroutine dist_to_boundary(neut, dist, surf, other_cell)
subroutine dist_to_boundary(p, dist, surf, other_cell)
type(Neutron), pointer :: neut
type(Particle), pointer :: p
real(8), intent(out) :: dist
integer, intent(out) :: surf
integer, optional, intent(in) :: other_cell
@ -246,7 +250,7 @@ contains
type(Cell), pointer :: parent_p => null()
type(Surface), pointer :: surf_p => null()
integer :: i
integer :: n_boundaries, n1, n2
integer :: n_surfaces, n1, n2
integer, allocatable :: expression(:)
integer :: index_surf
integer :: current_surf
@ -261,36 +265,36 @@ contains
if (present(other_cell)) then
cell_p => cells(other_cell)
else
cell_p => cells(neut%cell)
cell_p => cells(p%cell)
end if
current_surf = neut%surface
current_surf = p%surface
! determine number of surfaces to check
n1 = cell_p % n_items
n1 = cell_p % n_surfaces
n2 = 0
if (cell_p % parent > 0) then
parent_p => cells(cell_p % parent)
n2 = parent_p % n_items
n2 = parent_p % n_surfaces
end if
n_boundaries = n1 + n2
n_surfaces = n1 + n2
! allocate space and assign expression
allocate(expression(n_boundaries))
expression(1:n1) = cell_p%boundary_list
allocate(expression(n_surfaces))
expression(1:n1) = cell_p % surfaces
if (cell_p % parent > 0) then
expression(n1+1:n1+n2) = parent_p % boundary_list
expression(n1+1:n1+n2) = parent_p % surfaces
end if
! loop over all surfaces
dist = INFINITY
do i = 1, n_boundaries
x = neut%xyz(1)
y = neut%xyz(2)
z = neut%xyz(3)
u = neut%uvw(1)
v = neut%uvw(2)
w = neut%uvw(3)
do i = 1, n_surfaces
x = p % xyz(1)
y = p % xyz(2)
z = p % xyz(3)
u = p % uvw(1)
v = p % uvw(2)
w = p % uvw(3)
! check for coincident surface
index_surf = expression(i)
@ -306,7 +310,7 @@ contains
if (u == ZERO) then
d = INFINITY
else
x0 = surf_p%coeffs(1)
x0 = surf_p % coeffs(1)
d = (x0 - x)/u
if (d < ZERO) d = INFINITY
end if
@ -315,7 +319,7 @@ contains
if (v == ZERO) then
d = INFINITY
else
y0 = surf_p%coeffs(1)
y0 = surf_p % coeffs(1)
d = (y0 - y)/v
if (d < ZERO) d = INFINITY
end if
@ -324,16 +328,16 @@ contains
if (w == ZERO) then
d = INFINITY
else
z0 = surf_p%coeffs(1)
z0 = surf_p % coeffs(1)
d = (z0 - z)/w
if (d < ZERO) d = INFINITY
end if
case (SURF_PLANE)
A = surf_p%coeffs(1)
B = surf_p%coeffs(2)
C = surf_p%coeffs(3)
D = surf_p%coeffs(4)
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 == ZERO) then
@ -348,9 +352,9 @@ contains
if (a == ZERO) then
d = INFINITY
else
y0 = surf_p%coeffs(1)
z0 = surf_p%coeffs(2)
r = surf_p%coeffs(3)
y0 = surf_p % coeffs(1)
z0 = surf_p % coeffs(2)
r = surf_p % coeffs(3)
y = y - y0
z = z - z0
@ -383,9 +387,9 @@ contains
if (a == ZERO) then
d = INFINITY
else
x0 = surf_p%coeffs(1)
z0 = surf_p%coeffs(2)
r = surf_p%coeffs(3)
x0 = surf_p % coeffs(1)
z0 = surf_p % coeffs(2)
r = surf_p % coeffs(3)
x = x - x0
z = z - z0
@ -418,9 +422,9 @@ contains
if (a == ZERO) then
d = INFINITY
else
x0 = surf_p%coeffs(1)
y0 = surf_p%coeffs(2)
r = surf_p%coeffs(3)
x0 = surf_p % coeffs(1)
y0 = surf_p % coeffs(2)
r = surf_p % coeffs(3)
x = x - x0
y = y - y0
@ -454,10 +458,10 @@ contains
end if
case (SURF_SPHERE)
x0 = surf_p%coeffs(1)
y0 = surf_p%coeffs(2)
z0 = surf_p%coeffs(3)
r = surf_p%coeffs(4)
x0 = surf_p % coeffs(1)
y0 = surf_p % coeffs(2)
z0 = surf_p % coeffs(3)
r = surf_p % coeffs(4)
x = x - x0
y = y - y0
@ -531,56 +535,56 @@ contains
y = xyz(2)
z = xyz(3)
select case (surf%type)
select case (surf % type)
case (SURF_PX)
x0 = surf%coeffs(1)
x0 = surf % coeffs(1)
func = x - x0
case (SURF_PY)
y0 = surf%coeffs(1)
y0 = surf % coeffs(1)
func = y - y0
case (SURF_PZ)
z0 = surf%coeffs(1)
z0 = surf % coeffs(1)
func = z - z0
case (SURF_PLANE)
A = surf%coeffs(1)
B = surf%coeffs(2)
C = surf%coeffs(3)
D = surf%coeffs(4)
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)
y0 = surf%coeffs(1)
z0 = surf%coeffs(2)
r = surf%coeffs(3)
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)
x0 = surf%coeffs(1)
z0 = surf%coeffs(2)
r = surf%coeffs(3)
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)
x0 = surf%coeffs(1)
y0 = surf%coeffs(2)
r = surf%coeffs(3)
x0 = surf % coeffs(1)
y0 = surf % coeffs(2)
r = surf % coeffs(3)
func = (x-x0)**2 + (y-y0)**2 - r**2
case (SURF_SPHERE)
x0 = surf%coeffs(1)
y0 = surf%coeffs(2)
z0 = surf%coeffs(3)
r = surf%coeffs(4)
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)
y0 = surf%coeffs(1)
z0 = surf%coeffs(2)
y1 = surf%coeffs(3)
z1 = surf%coeffs(4)
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
s = SENSE_NEGATIVE
else
@ -589,10 +593,10 @@ contains
return
case (SURF_BOX_Y)
x0 = surf%coeffs(1)
z0 = surf%coeffs(2)
x1 = surf%coeffs(3)
z1 = surf%coeffs(4)
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
s = SENSE_NEGATIVE
else
@ -601,10 +605,10 @@ contains
return
case (SURF_BOX_Z)
x0 = surf%coeffs(1)
y0 = surf%coeffs(2)
x1 = surf%coeffs(3)
y1 = surf%coeffs(4)
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
s = SENSE_NEGATIVE
else
@ -613,12 +617,12 @@ contains
return
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)
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
s = SENSE_NEGATIVE
@ -672,8 +676,8 @@ contains
c => cells(i)
! loop over each surface specification
do j = 1, c%n_items
index = c%boundary_list(j)
do j = 1, c % n_surfaces
index = c % surfaces(j)
positive = (index > 0)
index = abs(index)
if (positive) then
@ -703,8 +707,8 @@ contains
c => cells(i)
! loop over each surface specification
do j = 1, c%n_items
index = c%boundary_list(j)
do j = 1, c % n_surfaces
index = c % surfaces(j)
positive = (index > 0)
index = abs(index)

10
src/global.f90
View file

@ -54,8 +54,8 @@ module global
type(ExtSource), target :: external_source
! Source and fission bank
type(Neutron), allocatable, target :: source_bank(:)
type(Bank), allocatable, target :: fission_bank(:)
type(Particle), allocatable, target :: source_bank(:)
type(Bank), allocatable, target :: fission_bank(:)
! Paths to input file, cross section data, etc
character(100) :: &
@ -142,9 +142,9 @@ module global
! Particle type
integer, parameter :: &
& NEUTRON_ = 1, &
& PHOTON_ = 2, &
& ELECTRON_ = 3
& NEUTRON = 1, &
& PHOTON = 2, &
& ELECTRON = 3
! Integer code for read error -- better hope this number is never
! used in an input file!

8
src/main.f90
View file

@ -90,7 +90,7 @@ contains
integer :: i, j
integer :: i_cycle
integer :: i_particle
type(Neutron), pointer :: particle => null()
type(Particle), pointer :: p => null()
CYCLE_LOOP: do i_cycle = 1, n_cycles
@ -99,8 +99,8 @@ contains
HISTORY_LOOP: do j = 1, n_particles
! grab source particle from bank
particle => get_source_particle()
if ( .not. associated(particle) ) then
p => get_source_particle()
if ( .not. associated(p) ) then
! no particles left in source bank
exit HISTORY_LOOP
end if
@ -110,7 +110,7 @@ contains
call RN_init_particle(int(i_particle,8))
! transport particle
call transport(particle)
call transport(p)
end do HISTORY_LOOP

6
src/output.f90
View file

@ -256,8 +256,8 @@ module output
end if
write(ou,*) ' Parent Cell = ' // int_to_str(c % parent)
string = ""
do i = 1, c % n_items
select case (c % boundary_List(i))
do i = 1, c % n_surfaces
select case (c % surfaces(i))
case (OP_LEFT_PAREN)
string = trim(string) // ' ('
case (OP_RIGHT_PAREN)
@ -267,7 +267,7 @@ module output
case (OP_DIFFERENCE)
string = trim(string) // ' !'
case default
string = trim(string) // ' ' // int_to_str(c % boundary_list(i))
string = trim(string) // ' ' // int_to_str(c % surfaces(i))
end select
end do
write(ou,*) ' Surface Specification:' // trim(string)

116
src/physics.f90
View file

@ -2,7 +2,7 @@ module physics
use global
use geometry, only: find_cell, dist_to_boundary, cross_boundary
use types, only: Neutron
use types, only: Particle
use mcnp_random, only: rang
use output, only: error, message
use search, only: binary_search
@ -17,9 +17,9 @@ contains
! geometry.
!=====================================================================
subroutine transport(neut)
subroutine transport(p)
type(Neutron), pointer :: neut
type(Particle), pointer :: p
character(250) :: msg
integer :: i, surf
@ -31,46 +31,45 @@ contains
real(8) :: distance
real(8) :: Sigma ! total cross-section
real(8) :: f ! interpolation factor
real(8) :: tmp(3)
integer :: IE ! index on energy grid
type(Universe), pointer :: univ
if (neut%cell == 0) then
if (p % cell == 0) then
univ => universes(BASE_UNIVERSE)
call find_cell(univ, neut, found_cell)
call find_cell(univ, p, found_cell)
! if neutron couldn't be located, print error
! if particle couldn't be located, print error
if (.not. found_cell) then
write(msg, 100) "Could not locate cell for neutron at: ", neut%xyz
write(msg, 100) "Could not locate cell for particle at: ", p % xyz
100 format (A,3ES11.3)
call error(msg)
end if
end if
if (verbosity >= 10) then
msg = "=== Particle " // trim(int_to_str(neut%uid)) // " ==="
msg = "=== Particle " // trim(int_to_str(p % uid)) // " ==="
call message(msg, 10)
i = cells(neut%cell)%uid
i = cells(p % cell)%uid
msg = " Born in cell " // trim(int_to_str(i))
call message(msg, 10)
end if
! sample energy from Watt fission energy spectrum for U-235
neut%E = watt_spectrum(0.988_8, 2.249_8)
p % E = watt_spectrum(0.988_8, 2.249_8)
! find energy index, interpolation factor
call find_energy_index(neut)
IE = neut % IE
f = neut % interp
call find_energy_index(p)
IE = p % IE
f = p % interp
! Determine material total cross-section
Sigma = f*cMaterial%total_xs(IE) + (1-f)*cMaterial%total_xs(IE+1)
do while (neut%alive)
do while (p % alive)
! Find the distance to the nearest boundary
call dist_to_boundary(neut, d_to_boundary, surf)
call dist_to_boundary(p, d_to_boundary, surf)
! Sample a distance to collision
d_to_collision = -log(rang()) / 1.0 ! Sigma
@ -78,18 +77,18 @@ contains
! Select smaller of the two distances
distance = min(d_to_boundary, d_to_collision)
! Advance neutron
neut%xyz = neut%xyz + distance*neut%uvw
! Advance particle
p%xyz = p%xyz + distance * p%uvw
! Add pathlength tallies
if (d_to_collision > d_to_boundary) then
neut%surface = surf
neut%cell = 0
call cross_boundary(neut)
p % surface = surf
p % cell = 0
call cross_boundary(p)
else
! collision
call collision(neut)
call collision(p)
end if
end do
@ -101,16 +100,16 @@ contains
! the interpolation factor for a particle at a certain energy
!=====================================================================
subroutine find_energy_index(neut)
subroutine find_energy_index(p)
type(Neutron), pointer :: neut
type(Particle), pointer :: p
integer :: IE
real(8) :: E
real(8) :: interp
! copy particle's energy
E = neut%E
E = p % E
! if particle's energy is outside of energy grid range, set to
! first or last index. Otherwise, do a binary search through the
@ -129,8 +128,8 @@ contains
interp = (E - e_grid(IE))/(e_grid(IE+1) - e_grid(IE))
! set particle attributes
neut % IE = IE
neut % interp = interp
p % IE = IE
p % interp = interp
end subroutine find_energy_index
@ -138,9 +137,9 @@ contains
! COLLISION
!=====================================================================
subroutine collision(neut)
subroutine collision(p)
type(Neutron), pointer :: neut
type(Particle), pointer :: p
type(AceContinuous), pointer :: table
type(AceReaction), pointer :: rxn
@ -151,7 +150,7 @@ contains
integer :: IE
real(8) :: f, Sigma, total
real(8) :: density, density_i
real(8) :: p
real(8) :: prob
real(8), allocatable :: Sigma_t(:)
! tallies
@ -167,8 +166,8 @@ contains
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))
IE = table%grid_index(p % IE)
f = (p%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
@ -176,18 +175,18 @@ contains
! sample nuclide
r1 = rang()
p = 0.0_8
prob = 0.0_8
total = sum(Sigma_t)
do i = 1, n_isotopes
p = p + Sigma_t(i) / total
if (r1 < p) exit
prob = prob + Sigma_t(i) / total
if (r1 < prob) exit
end do
! Get table, total xs, interpolation factor
table => xs_continuous(cMaterial%table(i))
Sigma = Sigma_t(i)
IE = table%grid_index(neut%IE)
f = (neut%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE))
IE = table%grid_index(p % IE)
f = (p%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE))
density = cMaterial%atom_percent(i)*density
! free memory
@ -195,13 +194,14 @@ contains
! sample reaction channel
r1 = rang()*Sigma
p = 0.0_8
prob = 0.0_8
do i = 1, table%n_reaction
rxn => table%reactions(i)
if (rxn%MT >= 200) cycle
if (IE < rxn%IE) cycle
p = p + density * (f*rxn%sigma(IE-rxn%IE+1) + (1-f)*(rxn%sigma(IE-rxn%IE+2)))
if (r1 < p) exit
prob = prob + density * (f*rxn%sigma(IE-rxn%IE+1) &
& + (1-f)*(rxn%sigma(IE-rxn%IE+2)))
if (r1 < prob) exit
end do
if (verbosity >= 10) then
msg = " " // trim(reaction_name(rxn%MT)) // " with nuclide " // &
@ -212,11 +212,11 @@ contains
! call appropriate subroutine
select case (rxn%MT)
case (2)
call elastic_scatter(neut, table%awr)
call elastic_scatter(p, table%awr)
case (102)
call n_gamma(neut)
call n_gamma(p)
case default
call elastic_scatter(neut, table%awr)
call elastic_scatter(p, table%awr)
end select
end subroutine collision
@ -225,9 +225,9 @@ contains
! ELASTIC_SCATTER
!=====================================================================
subroutine elastic_scatter(neut, awr)
subroutine elastic_scatter(p, awr)
type(Neutron), pointer :: neut
type(Particle), pointer :: p
real(8), intent(in) :: awr
real(8) :: phi ! azimuthal angle
@ -239,11 +239,11 @@ contains
real(8) :: E
integer :: IE
vel = sqrt(neut%E)
vel = sqrt(p % E)
vx = vel*neut%uvw(1)
vy = vel*neut%uvw(2)
vz = vel*neut%uvw(3)
vx = vel*p % uvw(1)
vy = vel*p % uvw(2)
vz = vel*p % uvw(3)
vcx = vx/(awr + 1.0_8)
@ -277,13 +277,13 @@ contains
E = vx*vx + vy*vy + vz*vz
vel = sqrt(E)
neut%E = E
neut%uvw(1) = vx/vel
neut%uvw(2) = vy/vel
neut%uvw(3) = vz/vel
p % E = E
p % uvw(1) = vx/vel
p % uvw(2) = vy/vel
p % uvw(3) = vz/vel
! find energy index, interpolation factor
call find_energy_index(neut)
call find_energy_index(p)
end subroutine elastic_scatter
@ -299,16 +299,16 @@ contains
! N_GAMMA
!=====================================================================
subroutine n_gamma(neut)
subroutine n_gamma(p)
type(Neutron), pointer :: neut
type(Particle), pointer :: p
integer :: cell_num
character(250) :: msg
neut%alive = .false.
p % alive = .false.
if (verbosity >= 10) then
cell_num = cells(neut%cell)%uid
cell_num = cells(p % cell)%uid
msg = " Absorbed in cell " // trim(int_to_str(cell_num))
call message(msg, 10)
end if

4
src/source.f90
View file

@ -65,13 +65,13 @@ contains
function get_source_particle()
type(Neutron), pointer :: get_source_particle
type(Particle), pointer :: get_source_particle
! increment index
source_index = source_index + 1
! if at end of bank, return null pointer
if ( source_index > size(source_bank) ) then
if (source_index > size(source_bank)) then
get_source_particle => null()
return
end if

63
src/types.f90
View file

@ -7,14 +7,14 @@ module types
!=====================================================================
type Universe
integer :: uid
integer :: type
integer :: level
integer :: n_cells
integer, allocatable :: cells(:)
real(8) :: x0
real(8) :: y0
real(8) :: z0
integer :: uid ! Unique ID
integer :: type ! Type
integer :: level ! Level of universe (0=base)
integer :: n_cells ! # of cells within
integer, allocatable :: cells(:) ! List of cells within
real(8) :: x0 ! Translation in x-coordinate
real(8) :: y0 ! Translation in y-coordinate
real(8) :: z0 ! Translation in z-coordinate
end type Universe
!=====================================================================
@ -22,10 +22,10 @@ module types
!=====================================================================
type Lattice
integer :: uid
integer :: type
integer :: level
real(8) :: pitch
integer :: uid ! Unique ID
integer :: type ! Type of lattice (square, hex, etc)
integer :: level ! Level of lattice
real(8) :: pitch ! Lattice pitch in cm
end type Lattice
!=====================================================================
@ -34,12 +34,13 @@ module types
!=====================================================================
type Surface
integer :: uid
integer :: type
real(8), allocatable :: coeffs(:)
integer, allocatable :: neighbor_pos(:)
integer, allocatable :: neighbor_neg(:)
integer :: bc
integer :: uid ! Unique ID
integer :: type ! Type of surface
real(8), allocatable :: coeffs(:) ! Definition of surface
integer, allocatable :: &
& neighbor_pos(:), & ! List of cells on positive side
& neighbor_neg(:) ! List of cells on negative side
integer :: bc ! Boundary condition
end type Surface
!=====================================================================
@ -47,23 +48,27 @@ module types
!=====================================================================
type Cell
integer :: uid
integer :: type
integer :: universe ! universe # this cell is in
integer :: fill ! universe # filling this cell
integer :: parent ! cell within which this cell resides
integer :: material
integer :: n_items
integer, allocatable :: boundary_list(:)
integer :: uid ! Unique ID
integer :: type ! Type of cell (normal, universe, lattice)
integer :: universe ! universe # this cell is in
integer :: fill ! universe # filling this cell
integer :: parent ! cell within which this cell resides
integer :: material ! Material within cell (0 for universe)
integer :: n_surfaces ! Number of surfaces within
integer, allocatable :: &
& surfaces(:) ! List of surfaces bounding cell -- note
! that parentheses, union, etc operators
! will be listed here too
end type Cell
!=====================================================================
! NEUTRON describes the state of a neutron being transported through
! PARTICLE describes the state of a particle being transported through
! the geometry
!=====================================================================
type Neutron
type Particle
integer :: uid ! Unique ID
integer :: type ! Particle type (n, p, e, etc)
real(8) :: xyz(3) ! location
real(8) :: uvw(3) ! directional cosines
real(8) :: E ! energy
@ -74,7 +79,7 @@ module types
integer :: surface ! current surface
real(8) :: wgt ! particle weight
logical :: alive ! is particle alive?
end type Neutron
end type Particle
!=====================================================================
! BANK is used for storing fission sites in criticality