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

Rectangular lattices now implemented.

Browse source
This commit is contained in:
Paul Romano 2011-03-28 20:56:06 +00:00
parent ea280fa9df
commit 9386d6fb10
9 changed files with 915 additions and 510 deletions
Download patch file Download diff file Expand all files Collapse all files

16
ChangeLog
View file

@ -1,3 +1,19 @@
2011-03-28 Paul Romano <romano7@mit.edu>
* fileio.f90: Ability to read lattice entries.
* geometry.f90: Added lattice capabilities. Changed cross_boundary
to cross_surface. Added many comments. In many places, changed
testing of c%fill to c%type instead since fill could be universe
or lattice.
* Makefile: Added -traceback
* output.f90: Added print_lattice, other minor changes. (Fixed
indentation so looks like everything changed)
* physics.f90: Call cross_lattice if crossing lattice boundary
instead of surface.
* source.f90: Set values for xyz_local, lattice indices
* types.f90: Extending Lattice type, modified Particle type to
include local coordinates, current lattice, and lattice indices.
2011-03-25 Paul Romano <romano7@mit.edu>
* physics.f90: Resample outgoing energy if sampled energy is

2
src/Makefile
View file

@ -26,7 +26,7 @@ test_objects = $(modules:.f90=.o) $(test:.f90=.o)
#--------------------------------------------------------------------
# Compiler Options
F90 = ifort
F90FLAGS = -g -fpp
F90FLAGS = -g -fpp -traceback
MPI = /opt/mpich2/1.3.3-gcc
MPIF90 = $(MPI)/bin/mpif90
MPI_COMPILE_FLAGS = $(shell mpif90 -f90='' -compile_info)

135
src/fileio.f90
View file

@ -97,6 +97,7 @@ contains
call dict_create(ucount_dict)
call dict_create(universe_dict)
call dict_create(lattice_dict)
open(FILE=filename, UNIT=in, STATUS='old', &
& ACTION='read', IOSTAT=ioError)
@ -131,6 +132,15 @@ contains
n_surfaces = n_surfaces + 1
case ('material')
n_materials = n_materials + 1
case ('lattice')
n_lattices = n_lattices + 1
! For lattices, 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(2))
call dict_add_key(lattice_dict, universe_num, n_lattices)
end select
end do
@ -155,6 +165,7 @@ contains
! Allocate arrays for cells, surfaces, etc.
allocate(cells(n_cells))
allocate(universes(n_universes))
allocate(lattices(n_lattices))
allocate(surfaces(n_surfaces))
allocate(materials(n_materials))
@ -251,7 +262,7 @@ contains
call read_cell(index_cell, words, n)
case ('surface')
! Read data
! Read data from surface entry
index_surface = index_surface + 1
call read_surface(index_surface, words, n)
@ -261,6 +272,7 @@ contains
call read_lattice(index_lattice, words, n)
case ('material')
! Read data from material entry
index_material = index_material + 1
call read_material(index_material, words, n)
@ -329,7 +341,7 @@ contains
end do
! adjust material indices
if (c % fill == 0) then
if (c % type == CELL_NORMAL) then
index = dict_get_key(material_dict, c%material)
if (index == DICT_NULL) then
msg = "Could not find material " // trim(int_to_str(c%material)) // &
@ -354,9 +366,13 @@ contains
integer, intent(in) :: level ! level of universe
integer :: i ! index for cells in universe
integer :: x,y ! indices for lattice positions
integer :: index ! index in cells array
integer :: universe_num
type(Cell), pointer :: c => null()
type(Universe), pointer :: subuniverse => null()
type(Lattice), pointer :: lat => null()
type(DictionaryII), pointer :: dict => null()
! set level of the universe
univ % level = level
@ -369,10 +385,31 @@ contains
! if this cell is filled with another universe, recursively
! call this subroutine
if (c % fill > 0) then
if (c % type == CELL_FILL) then
subuniverse => universes(c % fill)
call build_universe(subuniverse, index, level + 1)
end if
! if this cell is filled by a lattice, need to build the
! universe for each unique lattice element
if (c % type == CELL_LATTICE) then
lat => lattices(c % fill)
call dict_create(dict)
do x = 1, lat % n_x
do y = 1, lat % n_y
universe_num = lat % element(x,y)
if (dict_has_key(dict, universe_num)) then
cycle
else
call dict_add_key(dict, universe_num, 0)
subuniverse => universes(universe_num)
call build_universe(subuniverse, index, level + 1)
end if
end do
end do
end if
end do
end subroutine build_universe
@ -415,7 +452,6 @@ contains
! Read cell material
if (trim(words(4)) == 'fill') then
c % type = CELL_FILL
c % material = 0
n_surfaces = n_words - 5
@ -425,8 +461,15 @@ contains
msg = "Invalid universe fill: " // words(5)
call error(msg)
end if
c % fill = dict_get_key(universe_dict, universe_num)
! check whether universe or lattice
if (dict_has_key(universe_dict, universe_num)) then
c % type = CELL_FILL
c % fill = dict_get_key(universe_dict, universe_num)
elseif (dict_has_key(lattice_dict, universe_num)) then
c % type = CELL_LATTICE
c % fill = dict_get_key(lattice_dict, universe_num)
end if
else
c % type = CELL_NORMAL
c % material = str_to_int(words(4))
@ -439,7 +482,7 @@ contains
end if
! Assign number of items
c%n_surfaces = n_surfaces
c % n_surfaces = n_surfaces
! Read list of surfaces
allocate(c%surfaces(n_surfaces))
@ -560,9 +603,83 @@ contains
subroutine read_lattice(index, words, n_words)
integer, intent(in) :: index
character(*), intent(in) :: words(n_words)
integer, intent(in) :: n_words
integer, intent(in) :: index ! index in lattices array
character(*), intent(in) :: words(n_words) ! words in lattice entry
integer, intent(in) :: n_words ! number of words
integer :: universe_num ! user-specified universe number
integer :: n_x ! number of lattice cells in x direction
integer :: n_y ! number of lattice cells in y direction
integer :: i,j ! loop indices for Lattice % universes
integer :: index_word ! index in words array
character(250) :: msg ! output/error/message
character(32) :: word ! single word
type(Lattice), pointer :: lat => null()
lat => lattices(index)
! Read lattice universe
universe_num = str_to_int(words(2))
if (universe_num == ERROR_CODE) then
msg = "Invalid universe: " // words(2)
call error(msg)
end if
lat % uid = universe_num
! Read lattice type
word = words(3)
call lower_case(word)
select case(trim(word))
case ('rect')
lat % type = LATTICE_RECT
case ('hex')
lat % type = LATTICE_HEX
case default
msg = "Invalid lattice type: " // words(3)
call error(msg)
end select
! Read number of lattice cells in each direction
n_x = str_to_int(words(4))
n_y = str_to_int(words(5))
if (n_x == ERROR_CODE) then
msg = "Invalid number of lattice cells in x-direction: " // words(4)
call error(msg)
elseif (n_y == ERROR_CODE) then
msg = "Invalid number of lattice cells in y-direction: " // words(5)
call error(msg)
end if
lat % n_x = n_x
lat % n_y = n_y
! Read lattice origin coordinates
lat % x0 = str_to_real(words(6))
lat % y0 = str_to_real(words(7))
! Read lattice cell widths
lat % width_x = str_to_real(words(8))
lat % width_y = str_to_real(words(9))
! Make sure enough lattice positions specified
if (n_words - 9 < n_x * n_y) then
print *, n_words, n_x, n_y
msg = "Not enough lattice positions specified."
call error(msg)
end if
! Read lattice positions
allocate(lat % element(n_x,n_y))
do j = 0, n_y - 1
do i = 1, n_x
index_word = 9 + j*n_x + i
universe_num = str_to_int(words(index_word))
if (universe_num == ERROR_CODE) then
msg = "Invalid universe number: " // words(index_word)
call error(msg)
end if
lat % element(i, n_y-j) = dict_get_key(universe_dict, universe_num)
end do
end do
end subroutine read_lattice

368
src/geometry.f90
View file

@ -19,15 +19,15 @@ contains
type(Particle), pointer :: p
logical :: in_cell
integer, allocatable :: expression(:)
integer :: specified_sense
integer :: actual_sense
integer :: n_surfaces
integer :: i, j
integer :: surf_num
integer :: current_surface
integer, allocatable :: expression(:) ! copy of surfaces list
integer :: specified_sense ! specified sense of surface in list
integer :: actual_sense ! sense of particle wrt surface
integer :: n_surfaces ! number of surfaces in cell
integer :: i ! index of surfaces in cell
integer :: surf_num ! index in surfaces array (with sign)
integer :: current_surface ! current surface of particle (with sign)
character(250) :: msg ! output/error message
type(Surface), pointer :: surf => null()
character(250) :: msg
current_surface = p%surface
@ -59,7 +59,7 @@ contains
! Compare sense of point to specified sense
specified_sense = sign(1,expression(i))
call sense(surf, p%xyz, actual_sense)
actual_sense = sense(surf, p%xyz_local)
if (actual_sense == specified_sense) then
expression(i) = 1
else
@ -95,7 +95,9 @@ contains
character(250) :: msg ! error message
integer :: i ! index over cells
integer :: x, y
type(Cell), pointer :: c ! pointer to cell
type(Lattice), pointer :: lat ! pointer to lattice
type(Universe), pointer :: lower_univ ! if particle is in lower
! universe, use this pointer
! to call recursively
@ -107,9 +109,20 @@ contains
c => cells(univ % cells(i))
if (cell_contains(c, p)) then
! If this cell contains a universe of lattice, search for
! If this cell contains a universe or lattice, search for
! the particle in that universe/lattice
if (c % fill > 0) then
if (c % type == CELL_NORMAL) then
! set current pointers
found = .true.
cCell => c
cMaterial => materials(cCell%material)
cUniverse => univ
! set particle attributes
p%cell = univ % cells(i)
p%universe = dict_get_key(universe_dict, univ % uid)
exit
elseif (c % type == CELL_FILL) then
lower_univ => universes(c % fill)
call find_cell(lower_univ, p, found)
if (found) then
@ -118,17 +131,34 @@ contains
msg = "Could not locate particle in universe: "
call error(msg)
end if
else
! set current pointers
found = .true.
cCell => c
cMaterial => materials(cCell%material)
cUniverse => univ
elseif (c % type == CELL_LATTICE) then
! Set current lattice
lat => lattices(c % fill)
cLattice => lat
p % lattice = c % fill
! determine universe based on lattice position
x = ceiling((p%xyz(1) - lat%x0)/lat%width_x)
y = ceiling((p%xyz(2) - lat%y0)/lat%width_y)
lower_univ => universes(lat % element(x,y))
! adjust local position of particle
p%xyz_local(1) = p%xyz(1) - (lat%x0 + (x-0.5)*lat%width_x)
p%xyz_local(2) = p%xyz(2) - (lat%y0 + (y-0.5)*lat%width_y)
p%xyz_local(3) = p%xyz(3)
! set particle attributes
p%cell = dict_get_key(cell_dict, c % uid)
p%universe = dict_get_key(universe_dict, univ % uid)
exit
! set particle lattice indices
p % index_x = x
p % index_y = y
call find_cell(lower_univ, p, found)
if (found) then
exit
else
msg = "Could not locate particle in lattice: " &
& // int_to_str(lat % uid)
call error(msg)
end if
end if
end if
@ -137,19 +167,19 @@ contains
end subroutine find_cell
!=====================================================================
! CROSS_BOUNDARY moves a particle into a new cell
! CROSS_SURFACE moves a particle into a new cell
!=====================================================================
subroutine cross_boundary(p)
subroutine cross_surface(p)
type(Particle), pointer :: p
type(Surface), pointer :: surf
type(Cell), pointer :: c
integer :: i
integer :: index_cell
character(250) :: msg
logical :: found
integer :: i ! index of neighbors
integer :: index_cell ! index in cells array
logical :: found ! particle found in universe?
character(250) :: msg ! output/error message?
type(Surface), pointer :: surf
type(Cell), pointer :: c
type(Universe), pointer :: lower_univ => null()
surf => surfaces(abs(p%surface))
@ -175,7 +205,7 @@ contains
index_cell = surf%neighbor_pos(i)
c => cells(index_cell)
if (cell_contains(c, p)) then
if (c % fill > 0) then
if (c % type == CELL_FILL) then
lower_univ => universes(c % fill)
call find_cell(lower_univ, p, found)
if (.not. found) then
@ -198,7 +228,7 @@ contains
index_cell = surf%neighbor_neg(i)
c => cells(index_cell)
if (cell_contains(c, p)) then
if (c % fill > 0) then
if (c % type == CELL_FILL) then
lower_univ => universes(c % fill)
call find_cell(lower_univ, p, found)
if (.not. found) then
@ -233,44 +263,158 @@ contains
& ", it could not be located in any cell and it did not leak."
call error(msg)
end subroutine cross_boundary
end subroutine cross_surface
!=====================================================================
! DIST_TO_BOUNDARY calculates the distance to the nearest boundary of
! 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.
! CROSS_LATTICE moves a particle into a new lattice element
!=====================================================================
subroutine dist_to_boundary(p, dist, surf, other_cell)
subroutine cross_lattice(p)
type(Particle), pointer :: p
real(8), intent(out) :: dist
integer, intent(out) :: surf
integer, optional, intent(in) :: other_cell
type(Particle), pointer :: p
integer :: i_x ! x index in lattice
integer :: i_y ! y index in lattice
real(8) :: d_left ! distance to left side
real(8) :: d_right ! distance to right side
real(8) :: d_bottom ! distance to bottom side
real(8) :: d_top ! distance to top side
real(8) :: dist ! shortest distance
real(8) :: x ! x coordinate in local lattice element
real(8) :: y ! y coordinate in local lattice element
real(8) :: z ! z coordinate in local lattice element
real(8) :: u ! cosine of angle with x axis
real(8) :: v ! cosine of angle with y axis
real(8) :: x0 ! half the width of lattice element
real(8) :: y0 ! half the height of lattice element
logical :: found ! particle found in cell?
character(250) :: msg ! output/error message
type(Lattice), pointer :: lat
type(Universe), pointer :: univ
if (verbosity >= 10) then
msg = " Crossing lattice"
call message(msg, 10)
end if
lat => lattices(p % lattice)
u = p % uvw(1)
v = p % uvw(2)
if (lat % type == LATTICE_RECT) then
x = p % xyz_local(1)
y = p % xyz_local(2)
z = p % xyz_local(3)
x0 = lat % width_x * 0.5
y0 = lat % width_y * 0.5
dist = INFINITY
! left and right sides
if (u > 0) then
d_left = INFINITY
d_right = (x0 - x)/u
else
d_left = -(x0 + x)/u
d_right = INFINITY
end if
! top and bottom sides
if (v > 0) then
d_bottom = INFINITY
d_top = (y0 - y)/v
else
d_bottom = -(y0 + y)/v
d_top = INFINITY
end if
dist = min(d_left, d_right, d_top, d_bottom)
if (dist == d_left) then
! Move particle to left element
p % index_x = p % index_x - 1
p % xyz_local(1) = x0
elseif (dist == d_right) then
! Move particle to right element
p % index_x = p % index_x + 1
p % xyz_local(1) = -x0
elseif (dist == d_bottom) then
! Move particle to bottom element
p % index_y = p % index_y - 1
p % xyz_local(2) = y0
elseif (dist == d_top) then
! Move particle to top element
p % index_y = p % index_y + 1
p % xyz_local(2) = -y0
end if
elseif (lat % type == LATTICE_HEX) then
! TODO: Add hex lattice support
end if
! Check to make sure still in lattice
i_x = p % index_x
i_y = p % index_y
if (i_x < 1 .or. i_x > lat % n_x) then
msg = "Reached edge of lattice."
call error(msg)
elseif (i_y < 1 .or. i_y > lat % n_y) then
msg = "Reached edge of lattice."
call error(msg)
end if
! Find universe for next lattice element
univ => universes(lat % element(i_x,i_y))
! Find cell in next lattice element
call find_cell(univ, p, found)
if (.not. found) then
msg = "Could not locate particle in universe: "
call error(msg)
end if
end subroutine cross_lattice
!=====================================================================
! DIST_TO_BOUNDARY calculates the distance to the nearest boundary 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(p, dist, surf, in_lattice)
type(Particle), pointer :: p
real(8), intent(out) :: dist
integer, intent(out) :: surf
logical, intent(out) :: in_lattice
integer, allocatable :: expression(:) ! copy of surface list
integer :: i ! index for surface in cell
integer :: n_surfaces ! total number of surfaces to check
integer :: n1 ! number of surfaces in current cell
integer :: n2 ! number of surfaces in parent cell
integer :: index_surf ! index in surfaces array (with sign)
integer :: current_surf ! current surface
real(8) :: x,y,z ! particle coordinates
real(8) :: u,v,w ! particle directions
real(8) :: d ! evaluated distance
real(8) :: x0,y0,z0 ! coefficients for surface
real(8) :: r ! radius for quadratic surfaces
real(8) :: tmp ! dot product of surface normal with direction
real(8) :: a,b,c,k ! quadratic equation coefficients
real(8) :: quad ! discriminant of quadratic equation
logical :: on_surface ! is particle on surface?
character(250) :: msg ! output/error message
type(Cell), pointer :: cell_p => null()
type(Cell), pointer :: parent_p => null()
type(Surface), pointer :: surf_p => null()
integer :: i
integer :: n_surfaces, n1, n2
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
real(8) :: tmp
real(8) :: a,b,c,k
real(8) :: quad
character(250) :: msg
logical :: on_surface
type(LatticE), pointer :: lat => null()
if (present(other_cell)) then
cell_p => cells(other_cell)
else
cell_p => cells(p%cell)
end if
cell_p => cells(p%cell)
current_surf = p%surface
@ -290,16 +434,24 @@ contains
expression(n1+1:n1+n2) = parent_p % surfaces
end if
u = p % uvw(1)
v = p % uvw(2)
w = p % uvw(3)
! loop over all surfaces
dist = INFINITY
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)
if (i <= n1) then
! in local cell, so use xyz_local
x = p % xyz_local(1)
y = p % xyz_local(2)
z = p % xyz_local(3)
else
! in parent cell, so use xyz
x = p % xyz(1)
y = p % xyz(2)
z = p % xyz(3)
end if
! check for coincident surface
index_surf = expression(i)
if (index_surf == current_surf) then
@ -570,6 +722,45 @@ contains
end do
! Check lattice surfaces
in_lattice = .false.
if (p % lattice > 0) then
lat => lattices(p % lattice)
if (lat % type == LATTICE_RECT) then
x = p % xyz_local(1)
y = p % xyz_local(2)
z = p % xyz_local(3)
x0 = lat % width_x * 0.5
y0 = lat % width_y * 0.5
! left and right sides
if (u > 0) then
d = (x0 - x)/u
else
d = -(x + x0)/u
end if
if (d < dist) then
dist = d
in_lattice = .true.
end if
! top and bottom sides
if (v > 0) then
d = (y0 - y)/v
else
d = -(y + y0)/v
end if
if (d < dist) then
dist = d
in_lattice = .true.
end if
elseif (lat % type == LATTICE_HEX) then
! TODO: Add hex lattice support
end if
end if
! deallocate expression
deallocate(expression)
@ -581,17 +772,27 @@ contains
! a particular point is in.
!=====================================================================
subroutine sense(surf, xyz, s)
function sense(surf, xyz) result(s)
type(Surface), intent(in) :: surf
real(8), intent(in) :: xyz(3)
integer, intent(out) :: s
type(Surface), pointer :: surf ! surface
real(8), intent(in) :: xyz(3) ! coordinates of particle
integer :: s ! sense of particle
real(8) :: x,y,z
real(8) :: func
real(8) :: A,B,C,D,E,F,G,H,I,J
real(8) :: x0, y0, z0, r
real(8) :: x1, y1, z1
real(8) :: x,y,z ! coordinates of particle
real(8) :: func ! surface function evaluated at point
real(8) :: A ! coefficient on x**2 term in GQ
real(8) :: B ! coefficient on y**2 term in GQ
real(8) :: C ! coefficient on z**2 term in GQ
real(8) :: D ! coefficient on x*y term in GQ
real(8) :: E ! coefficient on y*z term in GQ
real(8) :: F ! coefficient on x*z term in GQ
real(8) :: G ! coefficient on x term in GQ
real(8) :: H ! coefficient on y term in GQ
real(8) :: I ! coefficient on z term in GQ
real(8) :: J ! coefficient on constant term in GQ
real(8) :: x0,y0,z0 ! coefficients for quadratic surfaces / box
real(8) :: r ! radius for quadratic surfaces
real(8) :: x1,y1,z1 ! upper-right corner of box
x = xyz(1)
y = xyz(2)
@ -715,7 +916,7 @@ contains
s = SENSE_NEGATIVE
end if
end subroutine sense
end function sense
!=====================================================================
! NEIGHBOR_LISTS builds a list of neighboring cells to each surface to
@ -724,16 +925,15 @@ contains
subroutine neighbor_lists()
integer :: i ! index in cells/surfaces array
integer :: j ! index of surface in cell
integer :: index ! index in count arrays
integer, allocatable :: count_positive(:) ! # of cells on positive side
integer, allocatable :: count_negative(:) ! # of cells on negative side
logical :: positive ! positive side specified in surface list
character(250) :: msg ! output/error message
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)

5
src/main.f90
View file

@ -5,7 +5,7 @@ program main
& normalize_ao, build_universe
use output, only: title, echo_input, message, warning, error, &
& print_summary, print_particle
use geometry, only: sense, cell_contains, neighbor_lists
use geometry, only: 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
@ -81,7 +81,10 @@ program main
call init_source()
! start problem
! surfaces(1)%bc = BC_VACUUM
! surfaces(2)%bc = BC_VACUUM
surfaces(1)%bc = BC_VACUUM
! surfaces(4)%bc = BC_VACUUM
call run_problem()
! deallocate arrays

821
src/output.f90
View file

@ -11,7 +11,7 @@ module output
integer :: ou = OUTPUT_UNIT
integer :: eu = ERROR_UNIT
contains
contains
!=====================================================================
! TITLE prints the main title banner as well as information about the
@ -19,131 +19,129 @@ module output
! run.
!=====================================================================
subroutine title()
subroutine title()
character(10) :: date
character(8) :: time
character(10) :: date
character(8) :: time
write(ou,*)
write(ou,*) ' .d88888b. 888b d888 .d8888b. '
write(ou,*) 'd88P" "Y88b 8888b d8888 d88P Y88b'
write(ou,*) '888 888 88888b.d88888 888 888'
write(ou,*) '888 888 88888b. .d88b. 88888b. 888Y88888P888 888 '
write(ou,*) '888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 '
write(ou,*) '888 888 888 888 88888888 888 888 888 Y8P 888 888 888'
write(ou,*) 'Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P'
write(ou,*) ' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P" '
write(ou,*) '____________888________________________________________________'
write(ou,*) ' 888 '
write(ou,*) ' 888 '
write(ou,*)
write(ou,*)
write(ou,*) ' .d88888b. 888b d888 .d8888b. '
write(ou,*) 'd88P" "Y88b 8888b d8888 d88P Y88b'
write(ou,*) '888 888 88888b.d88888 888 888'
write(ou,*) '888 888 88888b. .d88b. 88888b. 888Y88888P888 888 '
write(ou,*) '888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 '
write(ou,*) '888 888 888 888 88888888 888 888 888 Y8P 888 888 888'
write(ou,*) 'Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P'
write(ou,*) ' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P" '
write(ou,*) '____________888________________________________________________'
write(ou,*) ' 888 '
write(ou,*) ' 888 '
write(ou,*)
! Write version information
! write(ou,*) ' ______________________________________________________'
write(ou,*) ' Developed At: Massachusetts Institute of Technology'
write(ou,*) ' Lead Developer: Paul K. Romano'
write(ou,100) VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
100 format (6X,"Version:",9X,I1,".",I1,".",I1)
! Write version information
write(ou,*) ' Developed At: Massachusetts Institute of Technology'
write(ou,*) ' Lead Developer: Paul K. Romano'
write(ou,100) VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
100 format (6X,"Version:",9X,I1,".",I1,".",I1)
! Write the date and time
call get_today( date, time )
write(ou,101) trim(date), trim(time)
101 format (6X,"Date/Time:",7X,A,1X,A)
! write(ou,*) ' ______________________________________________________'
write(ou,*)
! Write the date and time
call get_today( date, time )
write(ou,101) trim(date), trim(time)
101 format (6X,"Date/Time:",7X,A,1X,A)
write(ou,*)
end subroutine title
end subroutine title
!=====================================================================
! ECHO_INPUT displays summary information about the problem about to
! be run after reading all the input.
!=====================================================================
subroutine echo_input()
subroutine echo_input()
character(32) :: string
character(32) :: string
! Display problem summary
write(ou,*) '============================================='
write(ou,*) '=> PROBLEM SUMMARY <='
write(ou,*) '============================================='
write(ou,*)
if ( problem_type == PROB_CRITICALITY ) then
write(ou,100) 'Problem type:', 'Criticality'
write(ou,100) 'Number of Cycles:', int_to_str(n_cycles)
write(ou,100) 'Number of Inactive Cycles:', int_to_str(n_inactive)
elseif ( problem_type == PROB_SOURCE ) then
write(ou,100) 'Problem type:', 'External Source'
end if
write(ou,100) 'Number of Particles:', int_to_str(n_particles)
write(ou,*)
! Display geometry summary
write(ou,*) '============================================='
write(ou,*) '=> GEOMETRY SUMMARY <='
write(ou,*) '============================================='
write(ou,*)
write(ou,100) 'Number of Cells:', int_to_str(n_cells)
write(ou,100) 'Number of Surfaces:', int_to_str(n_surfaces)
write(ou,100) 'Number of Materials:', int_to_str(n_materials)
write(ou,*)
! Display problem summary
write(ou,*) '============================================='
write(ou,*) '=> PROBLEM SUMMARY <='
write(ou,*) '============================================='
write(ou,*)
if ( problem_type == PROB_CRITICALITY ) then
write(ou,100) 'Problem type:', 'Criticality'
write(ou,100) 'Number of Cycles:', int_to_str(n_cycles)
write(ou,100) 'Number of Inactive Cycles:', int_to_str(n_inactive)
elseif ( problem_type == PROB_SOURCE ) then
write(ou,100) 'Problem type:', 'External Source'
end if
write(ou,100) 'Number of Particles:', int_to_str(n_particles)
write(ou,*)
! Display geometry summary
write(ou,*) '============================================='
write(ou,*) '=> GEOMETRY SUMMARY <='
write(ou,*) '============================================='
write(ou,*)
write(ou,100) 'Number of Cells:', int_to_str(n_cells)
write(ou,100) 'Number of Surfaces:', int_to_str(n_surfaces)
write(ou,100) 'Number of Materials:', int_to_str(n_materials)
write(ou,*)
! Format descriptor for columns
100 format (1X,A,T35,A)
! Format descriptor for columns
100 format (1X,A,T35,A)
end subroutine echo_input
end subroutine echo_input
!=====================================================================
! MESSAGE displays an informational message to the log file and the
! standard output stream.
!=====================================================================
subroutine message( msg, level )
subroutine message( msg, level )
character(*), intent(in) :: msg
integer, intent(in) :: level
character(*), intent(in) :: msg
integer, intent(in) :: level
integer :: n_lines
integer :: i
integer :: n_lines
integer :: i
! Only allow master to print to screen
if (.not. master) return
! Only allow master to print to screen
if (.not. master) return
if ( level <= verbosity ) then
n_lines = (len_trim(msg)-1)/79 + 1
do i = 1, n_lines
write(ou, fmt='(1X,A79)') msg(79*(i-1)+1:79*i)
end do
end if
if ( level <= verbosity ) then
n_lines = (len_trim(msg)-1)/79 + 1
do i = 1, n_lines
write(ou, fmt='(1X,A79)') msg(79*(i-1)+1:79*i)
end do
end if
end subroutine message
end subroutine message
!=====================================================================
! WARNING issues a warning to the user in the log file and the
! standard output stream.
!=====================================================================
subroutine warning( msg )
subroutine warning( msg )
character(*), intent(in) :: msg
character(*), intent(in) :: msg
integer :: n_lines
integer :: i
integer :: n_lines
integer :: i
! Only allow master to print to screen
if (.not. master) return
! Only allow master to print to screen
if (.not. master) return
write(ou, fmt='(1X,A9)', advance='no') 'WARNING: '
write(ou, fmt='(1X,A9)', advance='no') 'WARNING: '
n_lines = (len_trim(msg)-1)/70 + 1
do i = 1, n_lines
if ( i == 1 ) then
write(ou, fmt='(A70)') msg(70*(i-1)+1:70*i)
else
write(ou, fmt='(10X,A70)') msg(70*(i-1)+1:70*i)
end if
end do
end subroutine warning
n_lines = (len_trim(msg)-1)/70 + 1
do i = 1, n_lines
if ( i == 1 ) then
write(ou, fmt='(A70)') msg(70*(i-1)+1:70*i)
else
write(ou, fmt='(10X,A70)') msg(70*(i-1)+1:70*i)
end if
end do
end subroutine warning
!=====================================================================
! ERROR alerts the user that an error has been encountered and
@ -151,336 +149,369 @@ module output
! considered 'fatal' and hence the program is aborted.
!=====================================================================
subroutine error( msg )
subroutine error( msg )
character(*), intent(in) :: msg
character(*), intent(in) :: msg
integer :: n_lines
integer :: i
integer :: n_lines
integer :: i
! Only allow master to print to screen
if (master) then
write(eu, fmt='(1X,A7)', advance='no') 'ERROR: '
! Only allow master to print to screen
if (master) then
write(eu, fmt='(1X,A7)', advance='no') 'ERROR: '
n_lines = (len_trim(msg)-1)/72 + 1
do i = 1, n_lines
if ( i == 1 ) then
write(eu, fmt='(A72)') msg(72*(i-1)+1:72*i)
else
write(eu, fmt='(7X,A72)') msg(72*(i-1)+1:72*i)
end if
end do
write(eu,*)
end if
n_lines = (len_trim(msg)-1)/72 + 1
do i = 1, n_lines
if ( i == 1 ) then
write(eu, fmt='(A72)') msg(72*(i-1)+1:72*i)
else
write(eu, fmt='(7X,A72)') msg(72*(i-1)+1:72*i)
end if
end do
write(eu,*)
end if
! All processors abort
call free_memory()
end subroutine error
! All processors abort
call free_memory()
end subroutine error
!=====================================================================
! GET_TODAY determines the date and time at which the program began
! execution and returns it in a readable format
!=====================================================================
subroutine get_today(today_date, today_time)
subroutine get_today(today_date, today_time)
character(10), intent(out) :: today_date
character(8), intent(out) :: today_time
character(10), intent(out) :: today_date
character(8), intent(out) :: today_time
character(8) :: date
character(10) :: time
character(5) :: zone
integer :: val(8)
character(8) :: date
character(10) :: time
character(5) :: zone
integer :: val(8)
call date_and_time(date, time, zone, val)
! val(1) = year (YYYY)
! val(2) = month (MM)
! val(3) = day (DD)
! val(4) = timezone
! val(5) = hours (HH)
! val(6) = minutes (MM)
! val(7) = seconds (SS)
! val(8) = milliseconds
if ( val(2) < 10 ) then
if ( val(3) < 10 ) then
today_date = date(6:6) // "/" // date(8:8) // "/" // date(1:4)
else
today_date = date(6:6) // "/" // date(7:8) // "/" // date(1:4)
end if
else
if ( val(3) < 10 ) then
today_date = date(5:6) // "/" // date(8:8) // "/" // date(1:4)
else
today_date = date(5:6) // "/" // date(7:8) // "/" // date(1:4)
end if
end if
today_time = time(1:2) // ":" // time(3:4) // ":" // time(5:6)
call date_and_time(date, time, zone, val)
! val(1) = year (YYYY)
! val(2) = month (MM)
! val(3) = day (DD)
! val(4) = timezone
! val(5) = hours (HH)
! val(6) = minutes (MM)
! val(7) = seconds (SS)
! val(8) = milliseconds
end subroutine get_today
if ( val(2) < 10 ) then
if ( val(3) < 10 ) then
today_date = date(6:6) // "/" // date(8:8) // "/" // date(1:4)
else
today_date = date(6:6) // "/" // date(7:8) // "/" // date(1:4)
end if
else
if ( val(3) < 10 ) then
today_date = date(5:6) // "/" // date(8:8) // "/" // date(1:4)
else
today_date = date(5:6) // "/" // date(7:8) // "/" // date(1:4)
end if
end if
today_time = time(1:2) // ":" // time(3:4) // ":" // time(5:6)
end subroutine get_today
!=====================================================================
! PRINT_PARTICLE displays the attributes of a particle
!=====================================================================
subroutine print_particle(p)
subroutine print_particle(p)
type(Particle), pointer :: p
type(Particle), pointer :: p
integer :: i
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Universe), pointer :: u => null()
character(250) :: string
integer :: i
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Universe), pointer :: u => null()
character(250) :: string
select case (p % type)
case (NEUTRON)
write(ou,*) 'Neutron ' // int_to_str(p % uid)
case (PHOTON)
write(ou,*) 'Photon ' // int_to_str(p % uid)
case (ELECTRON)
write(ou,*) 'Electron ' // int_to_str(p % uid)
case default
write(ou,*) 'Unknown Particle ' // int_to_str(p % uid)
end select
write(ou,100) 'x = ', p % xyz(1)
write(ou,100) 'y = ', p % xyz(2)
write(ou,100) 'z = ', p % xyz(3)
write(ou,100) 'u = ', p % uvw(1)
write(ou,100) 'v = ', p % uvw(2)
write(ou,100) 'w = ', p % uvw(3)
write(ou,100) 'Weight = ', p % wgt
write(ou,100) 'Energy = ', p % E
write(ou,*) ' IE = ' // int_to_str(p % IE)
write(ou,100) 'Interpolation factor = ', p % interp
if (p % cell > 0) then
c => cells(p % cell)
write(ou,*) ' Cell = ' // int_to_str(c % uid)
else
write(ou,*) ' Cell not determined'
end if
select case (p % type)
case (NEUTRON)
write(ou,*) 'Neutron ' // int_to_str(p % uid)
case (PHOTON)
write(ou,*) 'Photon ' // int_to_str(p % uid)
case (ELECTRON)
write(ou,*) 'Electron ' // int_to_str(p % uid)
case default
write(ou,*) 'Unknown Particle ' // int_to_str(p % uid)
end select
write(ou,100) 'x = ', p % xyz(1)
write(ou,100) 'y = ', p % xyz(2)
write(ou,100) 'z = ', p % xyz(3)
write(ou,100) 'x local = ', p % xyz_local(1)
write(ou,100) 'y local = ', p % xyz_local(2)
write(ou,100) 'z local = ', p % xyz_local(3)
write(ou,100) 'u = ', p % uvw(1)
write(ou,100) 'v = ', p % uvw(2)
write(ou,100) 'w = ', p % uvw(3)
write(ou,100) 'Weight = ', p % wgt
write(ou,100) 'Energy = ', p % E
write(ou,*) ' x index = ' // int_to_str(p % index_x)
write(ou,*) ' y index = ' // int_to_str(p % index_y)
write(ou,*) ' IE = ' // int_to_str(p % IE)
write(ou,100) 'Interpolation factor = ', p % interp
if (p % surface > 0) then
s => surfaces(p % surface)
write(ou,*) ' Surface = ' // int_to_str(s % uid)
else
write(ou,*) ' Surface = None'
end if
if (p % cell > 0) then
c => cells(p % cell)
write(ou,*) ' Cell = ' // int_to_str(c % uid)
else
write(ou,*) ' Cell not determined'
end if
u => universes(p % universe)
write(ou,*) ' Universe = ' // int_to_str(u % uid)
write(ou,*)
if (p % surface > 0) then
s => surfaces(p % surface)
write(ou,*) ' Surface = ' // int_to_str(s % uid)
else
write(ou,*) ' Surface = None'
end if
! Format for a single real
100 format (5X,A,G10.3)
u => universes(p % universe)
write(ou,*) ' Universe = ' // int_to_str(u % uid)
write(ou,*)
nullify(c)
nullify(s)
nullify(u)
! Format for a single real
100 format (5X,A,G10.3)
end subroutine print_particle
nullify(c)
nullify(s)
nullify(u)
end subroutine print_particle
!=====================================================================
! PRINT_REACTION displays the attributes of a reaction
!=====================================================================
subroutine print_reaction(rxn)
subroutine print_reaction(rxn)
type(AceReaction), pointer :: rxn
type(AceReaction), pointer :: rxn
write(ou,*) 'Reaction ' // reaction_name(rxn % MT)
write(ou,*) ' MT = ' // int_to_str(rxn % MT)
write(ou,100) 'Q-value = ', rxn % Q_value
write(ou,*) ' TY = ' // int_to_str(rxn % TY)
write(ou,*) ' Starting index = ' // int_to_str(rxn % IE)
if (rxn % has_energy_dist) then
write(ou,*) ' Energy: Law ' // int_to_str(rxn % edist % law)
end if
write(ou,*)
write(ou,*) 'Reaction ' // reaction_name(rxn % MT)
write(ou,*) ' MT = ' // int_to_str(rxn % MT)
write(ou,100) 'Q-value = ', rxn % Q_value
write(ou,*) ' TY = ' // int_to_str(rxn % TY)
write(ou,*) ' Starting index = ' // int_to_str(rxn % IE)
if (rxn % has_energy_dist) then
write(ou,*) ' Energy: Law ' // int_to_str(rxn % edist % law)
end if
write(ou,*)
! Format for a single real
100 format (5X,A,G10.3)
! Format for a single real
100 format (5X,A,G10.3)
end subroutine print_reaction
end subroutine print_reaction
!=====================================================================
! PRINT_CELL displays the attributes of a cell
!=====================================================================
subroutine print_cell(c)
subroutine print_cell(c)
type(Cell), pointer :: c
type(Cell), pointer :: c
integer :: temp
integer :: i
type(Universe), pointer :: u => null()
type(Material), pointer :: m => null()
character(250) :: string
integer :: temp
integer :: i
type(Universe), pointer :: u => null()
type(Lattice), pointer :: l => null()
type(Material), pointer :: m => null()
character(250) :: string
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
write(ou,*) ' Parent Cell = ' // int_to_str(c % parent)
string = ""
do i = 1, c % n_surfaces
select case (c % surfaces(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 % surfaces(i))
end select
end do
write(ou,*) ' Surface Specification:' // trim(string)
write(ou,*)
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)
select case (c % type)
case (CELL_NORMAL)
write(ou,*) ' Fill = NONE'
case (CELL_FILL)
u => universes(c % fill)
write(ou,*) ' Fill = Universe ' // int_to_str(u % uid)
case (CELL_LATTICE)
l => lattices(c % fill)
write(ou,*) ' Fill = Lattice ' // int_to_str(l % uid)
end select
if (c % material == 0) then
write(ou,*) ' Material = NONE'
else
m => materials(c % material)
write(ou,*) ' Material = ' // int_to_str(m % uid)
end if
write(ou,*) ' Parent Cell = ' // int_to_str(c % parent)
string = ""
do i = 1, c % n_surfaces
select case (c % surfaces(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 % surfaces(i))
end select
end do
write(ou,*) ' Surface Specification:' // trim(string)
write(ou,*)
! nullify associated pointers
nullify(u)
nullify(m)
! nullify associated pointers
nullify(u)
nullify(m)
end subroutine print_cell
end subroutine print_cell
!=====================================================================
! PRINT_UNIVERSE displays the attributes of a universe
!=====================================================================
subroutine print_universe(univ)
subroutine print_universe(univ)
type(Universe), pointer :: univ
type(Universe), pointer :: univ
integer :: i
character(250) :: string
type(Cell), pointer :: c => null()
integer :: i
character(250) :: string
type(Cell), pointer :: c => null()
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,*)
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,*)
nullify(c)
nullify(c)
end subroutine print_universe
end subroutine print_universe
!=====================================================================
! PRINT_LATTICE displays the attributes of a lattice
!=====================================================================
subroutine print_lattice(lat)
type(Lattice), pointer :: lat
write(ou,*) 'Lattice ' // int_to_str(lat % uid)
write(ou,*) ' n_x = ' // int_to_str(lat % n_x)
write(ou,*) ' n_y = ' // int_to_str(lat % n_y)
write(ou,100) 'x0 = ', lat % x0
write(ou,100) 'y0 = ', lat % y0
write(ou,100) 'width_x = ', lat % width_x
write(ou,100) 'width_y = ', lat % width_y
write(ou,*)
! Format for a single real
100 format (5X,A,G10.3)
end subroutine print_lattice
!=====================================================================
! PRINT_SURFACE displays the attributes of a surface
!=====================================================================
subroutine print_surface(surf)
subroutine print_surface(surf)
type(Surface), pointer :: surf
type(Surface), pointer :: surf
integer :: i
character(80) :: string
integer :: i
character(80) :: string
write(ou,*) 'Surface ' // int_to_str(surf % uid)
select case (surf % type)
case (SURF_PX)
string = "X Plane"
case (SURF_PY)
string = "Y Plane"
case (SURF_PZ)
string = "Z Plane"
case (SURF_PLANE)
string = "Plane"
case (SURF_CYL_X)
string = "X Cylinder"
case (SURF_CYL_Y)
string = "Y Cylinder"
case (SURF_CYL_Z)
string = "Z Cylinder"
case (SURF_SPHERE)
string = "Sphere"
case (SURF_BOX_X)
case (SURF_BOX_Y)
case (SURF_BOX_Z)
case (SURF_BOX)
case (SURF_GQ)
string = "General Quadratic"
end select
write(ou,*) ' Type = ' // trim(string)
write(ou,*) ' Coefficients = ', surf % coeffs
write(ou,*) 'Surface ' // int_to_str(surf % uid)
select case (surf % type)
case (SURF_PX)
string = "X Plane"
case (SURF_PY)
string = "Y Plane"
case (SURF_PZ)
string = "Z Plane"
case (SURF_PLANE)
string = "Plane"
case (SURF_CYL_X)
string = "X Cylinder"
case (SURF_CYL_Y)
string = "Y Cylinder"
case (SURF_CYL_Z)
string = "Z Cylinder"
case (SURF_SPHERE)
string = "Sphere"
case (SURF_BOX_X)
case (SURF_BOX_Y)
case (SURF_BOX_Z)
case (SURF_BOX)
case (SURF_GQ)
string = "General Quadratic"
end select
write(ou,*) ' Type = ' // trim(string)
write(ou,*) ' Coefficients = ', surf % coeffs
string = ""
if (allocated(surf % neighbor_pos)) then
do i = 1, size(surf % neighbor_pos)
string = trim(string) // ' ' // int_to_str(surf % neighbor_pos(i))
end do
end if
write(ou,*) ' Positive Neighbors = ', trim(string)
string = ""
if (allocated(surf % neighbor_pos)) then
do i = 1, size(surf % neighbor_pos)
string = trim(string) // ' ' // int_to_str(surf % neighbor_pos(i))
end do
end if
write(ou,*) ' Positive Neighbors = ', trim(string)
string = ""
if (allocated(surf % neighbor_neg)) then
do i = 1, size(surf % neighbor_neg)
string = trim(string) // ' ' // int_to_str(surf % neighbor_neg(i))
end do
end if
write(ou,*) ' Negative Neighbors =', trim(string)
string = ""
if (allocated(surf % neighbor_neg)) then
do i = 1, size(surf % neighbor_neg)
string = trim(string) // ' ' // int_to_str(surf % neighbor_neg(i))
end do
end if
write(ou,*) ' Negative Neighbors =', trim(string)
write(ou,*)
end subroutine print_surface
end subroutine print_surface
!=====================================================================
! PRINT_MATERIAL displays the attributes of a material
!=====================================================================
subroutine print_material(mat)
subroutine print_material(mat)
type(Material), pointer :: mat
type(Material), pointer :: mat
integer :: i
integer :: n_lines
type(AceContinuous), pointer :: table
character(250) :: string
integer :: i
integer :: n_lines
type(AceContinuous), pointer :: table
character(250) :: string
write(ou,*) 'Material ' // int_to_str(mat % uid)
! Make string of all isotopes
string = ""
do i = 1, mat % n_isotopes
table => xs_continuous(mat % table(i))
string = trim(string) // ' ' // table % name
end do
! Print isotopes with word wrap
! TODO: Change this to generic word wrap subroutine?
n_lines = (len_trim(string)-1)/75 + 1
do i = 1, n_lines
if ( i == 1 ) then
write(ou, fmt='(5X,A75)') 'Isotopes =' // string(75*(i-1)+1:75*i)
else
write(ou, fmt='(5X,A75)') string(75*(i-1)+1:75*i)
end if
end do
write(ou,'(5X,A,G12.4,A)') 'Atom Density = ', mat % atom_density, &
& ' atom/b-cm'
write(ou,*)
write(ou,*) 'Material ' // int_to_str(mat % uid)
! Make string of all isotopes
string = ""
do i = 1, mat % n_isotopes
table => xs_continuous(mat % table(i))
string = trim(string) // ' ' // table % name
end do
! Print isotopes with word wrap
! TODO: Change this to generic word wrap subroutine?
n_lines = (len_trim(string)-1)/75 + 1
do i = 1, n_lines
if ( i == 1 ) then
write(ou, fmt='(5X,A75)') 'Isotopes =' // string(75*(i-1)+1:75*i)
else
write(ou, fmt='(5X,A75)') string(75*(i-1)+1:75*i)
end if
end do
write(ou,'(5X,A,G12.4,A)') 'Atom Density = ', mat % atom_density, &
& ' atom/b-cm'
write(ou,*)
nullify(table)
nullify(table)
end subroutine print_material
end subroutine print_material
!=====================================================================
! PRINT_SUMMARY displays the attributes of all cells, universes,
@ -488,59 +519,73 @@ module output
! debugging!
!=====================================================================
subroutine print_summary()
subroutine print_summary()
type(Surface), pointer :: s
type(Cell), pointer :: c
type(Universe), pointer :: u
type(Material), pointer :: m
integer :: i
type(Surface), pointer :: s => null()
type(Cell), pointer :: c => null()
type(Universe), pointer :: u => null()
type(Lattice), pointer :: l => null()
type(Material), pointer :: m => null()
integer :: i
! 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 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 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
write(ou,*) '============================================='
write(ou,*) '=> SURFACE SUMMARY <='
write(ou,*) '============================================='
write(ou,*)
do i = 1, n_surfaces
s => surfaces(i)
call print_surface(s)
end do
! print summary of lattices
if (n_lattices > 0) then
write(ou,*) '============================================='
write(ou,*) '=> LATTICE SUMMARY <='
write(ou,*) '============================================='
write(ou,*)
do i = 1, n_lattices
l => lattices(i)
call print_lattice(l)
end do
end if
! print summary of materials
write(ou,*) '============================================='
write(ou,*) '=> MATERIAL SUMMARY <='
write(ou,*) '============================================='
write(ou,*)
do i = 1, n_materials
m => materials(i)
call print_material(m)
end do
! print summary of surfaces
write(ou,*) '============================================='
write(ou,*) '=> SURFACE SUMMARY <='
write(ou,*) '============================================='
write(ou,*)
do i = 1, n_surfaces
s => surfaces(i)
call print_surface(s)
end do
nullify(s)
nullify(c)
nullify(u)
nullify(m)
! print summary of materials
write(ou,*) '============================================='
write(ou,*) '=> MATERIAL SUMMARY <='
write(ou,*) '============================================='
write(ou,*)
do i = 1, n_materials
m => materials(i)
call print_material(m)
end do
end subroutine print_summary
nullify(s)
nullify(c)
nullify(u)
nullify(l)
nullify(m)
end subroutine print_summary
end module output

16
src/physics.f90
View file

@ -1,7 +1,8 @@
module physics
use global
use geometry, only: find_cell, dist_to_boundary, cross_boundary
use geometry, only: find_cell, dist_to_boundary, cross_surface, &
& cross_lattice
use types, only: Particle
use mcnp_random, only: rang
use output, only: error, warning, message, print_particle
@ -29,6 +30,7 @@ contains
real(8) :: Sigma ! total cross-section
real(8) :: f ! interpolation factor
logical :: found_cell ! found cell which particle is in?
logical :: in_lattice ! is surface crossing in lattice?
character(250) :: msg ! output/error message
type(Universe), pointer :: univ
@ -67,7 +69,7 @@ contains
Sigma = (1-f)*cMaterial%total_xs(IE) + f*cMaterial%total_xs(IE+1)
! Find the distance to the nearest boundary
call dist_to_boundary(p, d_to_boundary, surf)
call dist_to_boundary(p, d_to_boundary, surf, in_lattice)
! Sample a distance to collision
d_to_collision = -log(rang()) / Sigma
@ -77,13 +79,19 @@ contains
! Advance particle
p%xyz = p%xyz + distance * p%uvw
p%xyz_local = p%xyz_local + distance * p%uvw
! Add pathlength tallies
if (d_to_collision > d_to_boundary) then
p % surface = surf
p % cell = 0
call cross_boundary(p)
if (in_lattice) then
p % surface = 0
call cross_lattice(p)
else
p % surface = surf
call cross_surface(p)
end if
else
! collision
p % surface = 0

13
src/source.f90
View file

@ -65,6 +65,7 @@ contains
r = (/ (rang(), k = 1,3) /)
p % uid = j
p % xyz = p_min + r*(p_max - p_min)
p % xyz_local = p % xyz
! angle
phi = TWO*PI*rang()
@ -78,8 +79,11 @@ contains
p % cell = 0
p % surface = 0
p % universe = 0
p % lattice = 0
p % wgt = ONE
p % alive = .true.
p % index_x = 0
p % index_y = 0
! sample energy from Watt fission energy spectrum for U-235
do
@ -160,10 +164,11 @@ contains
do i = 1, n_sites
j = index + i - 1
source_bank(j) % uid = temp_bank(i) % uid
source_bank(j) % xyz = temp_bank(i) % xyz
source_bank(j) % uvw = temp_bank(i) % uvw
source_bank(j) % E = temp_bank(i) % E
source_bank(j) % uid = temp_bank(i) % uid
source_bank(j) % xyz = temp_bank(i) % xyz
source_bank(j) % xyz_local = temp_bank(i) % xyz
source_bank(j) % uvw = temp_bank(i) % uvw
source_bank(j) % E = temp_bank(i) % E
! set defaults
source_bank(j) % type = NEUTRON

49
src/types.f90
View file

@ -3,7 +3,7 @@ module types
implicit none
!=====================================================================
! UNIVERSE
! UNIVERSE defines a geometry that fills all phase space
!=====================================================================
type Universe
@ -18,14 +18,21 @@ module types
end type Universe
!=====================================================================
! LATTICE
! LATTICE is an ordered array of elements (either rectangular,
! hexagonal, or triangular)
!=====================================================================
type Lattice
integer :: uid ! Unique ID
integer :: type ! Type of lattice (square, hex, etc)
integer :: level ! Level of lattice
real(8) :: pitch ! Lattice pitch in cm
integer :: uid ! Universe number for lattice
integer :: type ! Type of lattice (rectangular, hex, etc)
integer :: level ! Level of lattice
integer :: n_x ! number of lattice cells in x-direction
integer :: n_y ! number of lattice cells in y-direction
real(8) :: x0 ! x-coordinate of lattice origin
real(8) :: y0 ! y-coordinate of lattice origin
real(8) :: width_x ! width of lattice cell
real(8) :: width_y ! width of lattice cell
integer, allocatable :: element(:,:) ! specified universes
end type Lattice
!=====================================================================
@ -67,19 +74,23 @@ module types
!=====================================================================
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) :: wgt ! particle weight
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
logical :: alive ! is particle alive?
integer :: n_coll ! # of collisions
integer :: uid ! Unique ID
integer :: type ! Particle type (n, p, e, etc)
real(8) :: xyz(3) ! location
real(8) :: xyz_local(3) ! local location (after transformations)
real(8) :: uvw(3) ! directional cosines
real(8) :: wgt ! particle weight
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 :: lattice ! current lattice
integer :: surface ! current surface
integer :: index_x ! lattice index for x direction
integer :: index_y ! lattice index for y direction
logical :: alive ! is particle alive?
integer :: n_coll ! # of collisions
end type Particle
!=====================================================================