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

merged Paul's latest commits and re-compiled

Browse source
This commit is contained in:
Bryan Herman 2011-12-26 15:00:52 -08:00
commit a43db54057
31 changed files with 2850 additions and 2439 deletions
Download patch file Download diff file Expand all files Collapse all files

70
src/DEPENDENCIES
View file

@ -1,3 +1,18 @@
ace.o: ace_header.o
ace.o: constants.o
ace.o: datatypes.o
ace.o: datatypes_header.o
ace.o: endf.o
ace.o: error.o
ace.o: fission.o
ace.o: global.o
ace.o: material_header.o
ace.o: output.o
ace.o: string.o
ace_header.o: constants.o
ace_header.o: endf_header.o
cmfd_execute.o: cmfd_utils.o
cmfd_execute.o: global.o
cmfd_execute.o: mesh.o
@ -13,21 +28,15 @@ cmfd_utils.o: mesh_header.o
cmfd_utils.o: string.o
cmfd_utils.o: xml-fortran/templates/cmfd_t.o
cross_section.o: ace_header.o
cross_section.o: constants.o
cross_section.o: cross_section_header.o
cross_section.o: datatypes.o
cross_section.o: datatypes_header.o
cross_section.o: endf.o
cross_section.o: error.o
cross_section.o: fileio.o
cross_section.o: fission.o
cross_section.o: global.o
cross_section.o: global.o
cross_section.o: material_header.o
cross_section.o: output.o
cross_section.o: string.o
cross_section_header.o: constants.o
cross_section_header.o: endf_header.o
cross_section.o: particle_header.o
cross_section.o: random_lcg.o
cross_section.o: search.o
datatypes.o: datatypes_header.o
@ -44,12 +53,8 @@ energy_grid.o: datatypes_header.o
energy_grid.o: global.o
energy_grid.o: output.o
fileio.o: constants.o
fileio.o: global.o
fileio.o: string.o
fission.o: ace_header.o
fission.o: constants.o
fission.o: cross_section_header.o
fission.o: error.o
fission.o: global.o
fission.o: interpolation.o
@ -65,10 +70,10 @@ geometry.o: particle_header.o
geometry.o: string.o
geometry.o: tally.o
global.o: ace_header.o
global.o: bank_header.o
global.o: cmfd_header.o
global.o: constants.o
global.o: cross_section_header.o
global.o: datatypes_header.o
global.o: geometry_header.o
global.o: material_header.o
@ -78,8 +83,9 @@ global.o: source_header.o
global.o: tally_header.o
global.o: timing.o
initialize.o: ace.o
initialize.o: bank_header.o
initialize.o: constants.o
initialize.o: cross_section.o
initialize.o: datatypes.o
initialize.o: datatypes_header.o
initialize.o: energy_grid.o
@ -88,8 +94,6 @@ initialize.o: geometry.o
initialize.o: geometry_header.o
initialize.o: global.o
initialize.o: input_xml.o
initialize.o: logging.o
initialize.o: mpi_routines.o
initialize.o: output.o
initialize.o: random_lcg.o
initialize.o: source.o
@ -113,6 +117,14 @@ input_xml.o: xml-fortran/templates/materials_t.o
input_xml.o: xml-fortran/templates/settings_t.o
input_xml.o: xml-fortran/templates/tallies_t.o
intercycle.o: error.o
intercycle.o: global.o
intercycle.o: output.o
intercycle.o: particle_header.o
intercycle.o: random_lcg.o
intercycle.o: tally_header.o
intercycle.o: timing.o
interpolation.o: constants.o
interpolation.o: endf_header.o
interpolation.o: error.o
@ -127,7 +139,7 @@ main.o: cmfd_execute.o
main.o: constants.o
main.o: global.o
main.o: initialize.o
main.o: mpi_routines.o
main.o: intercycle.o
main.o: output.o
main.o: particle_header.o
main.o: physics.o
@ -140,15 +152,6 @@ main.o: timing.o
mesh.o: mesh_header.o
mpi_routines.o: constants.o
mpi_routines.o: error.o
mpi_routines.o: global.o
mpi_routines.o: output.o
mpi_routines.o: particle_header.o
mpi_routines.o: random_lcg.o
mpi_routines.o: tally_header.o
mpi_routines.o: timing.o
output.o: constants.o
output.o: datatypes.o
output.o: endf.o
@ -162,8 +165,9 @@ output.o: tally_header.o
particle_header.o: constants.o
particle_header.o: geometry_header.o
physics.o: ace_header.o
physics.o: constants.o
physics.o: cross_section_header.o
physics.o: cross_section.o
physics.o: endf.o
physics.o: error.o
physics.o: fission.o
@ -189,9 +193,9 @@ search.o: constants.o
search.o: error.o
search.o: global.o
source.o: ace_header.o
source.o: bank_header.o
source.o: constants.o
source.o: cross_section_header.o
source.o: error.o
source.o: global.o
source.o: output.o
@ -207,9 +211,9 @@ string.o: global.o
tally.o: constants.o
tally.o: error.o
tally.o: global.o
tally.o: intercycle.o
tally.o: mesh.o
tally.o: mesh_header.o
tally.o: mpi_routines.o
tally.o: output.o
tally.o: search.o
tally.o: string.o

18
src/Makefile
View file

@ -118,7 +118,7 @@ endif
ifeq ($(COMPILER),ibm)
F90 = xlf2003
F90FLAGS := -WF,-DNO_F2008
F90FLAGS := -WF,-DNO_F2008 -O2
# Debugging
ifeq ($(DEBUG),yes)
@ -166,11 +166,19 @@ ifeq ($(USE_MPI),yes)
endif
#===============================================================================
# Special options for ORNL Jaguar supercomputer
# Options for IBM Blue Gene/P ANL supercomputer
#===============================================================================
HOSTNAME = $(shell hostname)
ifneq (,$(findstring jaguar,$(HOSTNAME)))
ifeq ($(MACHINE),bluegene)
F90 = /bgsys/drivers/ppcfloor/comm/xl/bin/mpixlf2003
F90FLAGS = -WF,-DNO_F2008,-DMPI -O3
endif
#===============================================================================
# Options for Cray XK6 ORNL Jaguar supercomputer
#===============================================================================
ifeq ($(MACHINE),crayxk6)
F90 = ftn
F90FLAGS += -DMPI
endif
@ -181,7 +189,7 @@ endif
all: xml-fortran $(program)
xml-fortran:
cd xml-fortran; make F90=$(F90) F90FLAGS="$(F90FLAGS)"
cd xml-fortran; make MACHINE=$(MACHINE) F90=$(F90) F90FLAGS="$(F90FLAGS)"
cd xml-fortran/templates; make F90=$(F90) F90FLAGS="$(F90FLAGS)"
$(program): $(objects)
$(F90) $(objects) $(templates) $(xml_fort) -o $@ $(LDFLAGS)

7
src/OBJECTS
View file

@ -1,9 +1,10 @@
objects = \
ace.o \
ace_header.o \
bank_header.o \
cmfd_execute.o \
cmfd_header.o \
cmfd_utils.o \
cross_section_header.o \
cross_section.o \
datatypes.o \
datatypes_header.o \
@ -12,20 +13,18 @@ endf.o \
endf_header.o \
energy_grid.o \
error.o \
fileio.o \
fission.o \
geometry.o \
geometry_header.o \
global.o \
initialize.o \
intercycle.o \
interpolation.o \
input_xml.o \
logging.o \
main.o \
material_header.o \
mesh_header.o \
mesh.o \
mpi_routines.o \
output.o \
particle_header.o \
physics.o \

1314
src/ace.F90 Normal file
View file

File diff suppressed because it is too large Load diff

14
src/cross_section_header.F90 → src/ace_header.F90
View file

@ -1,4 +1,4 @@
module cross_section_header
module ace_header
use constants, only: MAX_FILE_LEN
use endf_header, only: Tab1
@ -51,11 +51,16 @@ module cross_section_header
end type Reaction
!===============================================================================
! URRDATA contains unresolved resonance data.
! URRDATA contains probability tables for the unresolved resonance range.
!===============================================================================
type UrrData
integer, allocatable :: params(:)
integer :: n_energy ! # of incident neutron energies
integer :: n_prob ! # of probabilities
integer :: interp ! inteprolation (2=lin-lin, 5=log-log)
integer :: inelastic_flag ! inelastic competition flag
integer :: absorption_flag ! other absorption flag
logical :: multiply_smooth ! multiply by smooth cross section?
real(8), allocatable :: energy(:)
real(8), allocatable :: prob(:,:,:)
end type UrrData
@ -108,6 +113,7 @@ module cross_section_header
! Unresolved resonance data
logical :: urr_present
integer :: urr_inelastic
type(UrrData), pointer :: urr_data => null()
! Reactions
@ -203,4 +209,4 @@ module cross_section_header
real(8) :: nu_fission ! macroscopic production xs
end type MaterialMacroXS
end module cross_section_header
end module ace_header

6
src/cmfd_utils.F90
View file

@ -395,9 +395,9 @@ contains
res = leakage + interactions - scattering - (ONE/keff)*fission
! write output
label = "MESH (" // trim(int_to_str(i)) // ". " // &
& trim(int_to_str(j)) // ", " // trim(int_to_str(k)) // &
& ") GROUP " // trim(int_to_str(g))
label = "MESH (" // trim(int4_to_str(i)) // ". " // &
& trim(int4_to_str(j)) // ", " // trim(int4_to_str(k)) // &
& ") GROUP " // trim(int4_to_str(g))
write(UNIT=UNIT_CMFD, FMT='(A,T35,A)') label, &
& trim(real_to_str(res))

112
src/constants.F90
View file

@ -22,7 +22,7 @@ module constants
real(8), parameter :: TINY_BIT = 1e-8_8
! User for precision in geometry
real(8), parameter :: FP_PRECISION = 1e-7_8
real(8), parameter :: FP_PRECISION = 1e-5_8
! Maximum number of collisions/crossings
integer, parameter :: MAX_EVENTS = 10000
@ -54,57 +54,64 @@ module constants
! GEOMETRY-RELATED CONSTANTS
! Boundary conditions
integer, parameter :: &
& BC_TRANSMIT = 0, & ! Transmission boundary condition (default)
& BC_VACUUM = 1, & ! Vacuum boundary condition
& BC_REFLECT = 2, & ! Reflecting boundary condition
& BC_PERIODIC = 3 ! Periodic boundary condition
integer, parameter :: &
BC_TRANSMIT = 0, & ! Transmission boundary condition (default)
BC_VACUUM = 1, & ! Vacuum boundary condition
BC_REFLECT = 2, & ! Reflecting boundary condition
BC_PERIODIC = 3 ! Periodic boundary condition
! Logical operators for cell definitions
integer, parameter :: &
& OP_LEFT_PAREN = huge(0), & ! Left parentheses
& OP_RIGHT_PAREN = huge(0) - 1, & ! Right parentheses
& OP_UNION = huge(0) - 2, & ! Union operator
& OP_DIFFERENCE = huge(0) - 3 ! Difference operator
integer, parameter :: &
OP_LEFT_PAREN = huge(0), & ! Left parentheses
OP_RIGHT_PAREN = huge(0) - 1, & ! Right parentheses
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
& CELL_VOID = -1
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
CELL_VOID = -1
! Lattice types
integer, parameter :: &
& LATTICE_RECT = 1, &
& LATTICE_HEX = 2
integer, parameter :: &
LATTICE_RECT = 1, &
LATTICE_HEX = 2
! Lattice boundary crossings
integer, parameter :: &
LATTICE_LEFT = 1, &
LATTICE_RIGHT = 2, &
LATTICE_BOTTOM = 3, &
LATTICE_TOP = 4
! Surface types
integer, parameter :: &
& SURF_PX = 1, & ! Plane parallel to x-plane
& SURF_PY = 2, & ! Plane parallel to y-plane
& SURF_PZ = 3, & ! Plane parallel to z-plane
& SURF_PLANE = 4, & ! Arbitrary plane
& SURF_CYL_X = 5, & ! Cylinder along x-axis
& SURF_CYL_Y = 6, & ! Cylinder along y-axis
& SURF_CYL_Z = 7, & ! Cylinder along z-axis
& SURF_SPHERE = 8, & ! Sphere
& SURF_BOX_X = 9, & ! Box extending infinitely in x-direction
& SURF_BOX_Y = 10, & ! Box extending infinitely in y-direction
& SURF_BOX_Z = 11, & ! Box extending infinitely in z-direction
& SURF_BOX = 12, & ! Rectangular prism
& SURF_GQ = 13 ! General quadratic surface
integer, parameter :: &
SURF_PX = 1, & ! Plane parallel to x-plane
SURF_PY = 2, & ! Plane parallel to y-plane
SURF_PZ = 3, & ! Plane parallel to z-plane
SURF_PLANE = 4, & ! Arbitrary plane
SURF_CYL_X = 5, & ! Cylinder along x-axis
SURF_CYL_Y = 6, & ! Cylinder along y-axis
SURF_CYL_Z = 7, & ! Cylinder along z-axis
SURF_SPHERE = 8, & ! Sphere
SURF_BOX_X = 9, & ! Box extending infinitely in x-direction
SURF_BOX_Y = 10, & ! Box extending infinitely in y-direction
SURF_BOX_Z = 11, & ! Box extending infinitely in z-direction
SURF_BOX = 12, & ! Rectangular prism
SURF_GQ = 13 ! General quadratic surface
! Surface senses
integer, parameter :: &
& SENSE_POSITIVE = 1, &
& SENSE_NEGATIVE = -1
integer, parameter :: &
SENSE_POSITIVE = 1, &
SENSE_NEGATIVE = -1
! ============================================================================
! CROSS SECTION RELATED CONSTANTS
! Interpolation flag
integer, parameter :: &
integer, parameter :: &
HISTOGRAM = 1, & ! y is constant in x
LINEAR_LINEAR = 2, & ! y is linear in x
LINEAR_LOG = 3, & ! y is linear in ln(x)
@ -119,19 +126,19 @@ module constants
! Angular distribution type
integer, parameter :: &
& ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution
& ANGLE_32_EQUI = 2, & ! 32 equiprobable bins
& ANGLE_TABULAR = 3 ! Tabular angular distribution
ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution
ANGLE_32_EQUI = 2, & ! 32 equiprobable bins
ANGLE_TABULAR = 3 ! Tabular angular distribution
! Secondary energy mode for S(a,b) inelastic scattering
integer, parameter :: &
& SAB_SECONDARY_EQUAL = 0, & ! Equally-likely outgoing energy bins
& SAB_SECONDARY_SKEWED = 1 ! Skewed outgoing energy bins
SAB_SECONDARY_EQUAL = 0, & ! Equally-likely outgoing energy bins
SAB_SECONDARY_SKEWED = 1 ! Skewed outgoing energy bins
! Elastic mode for S(a,b) elastic scattering
integer, parameter :: &
& SAB_ELASTIC_DISCRETE = 3, & ! Sample from discrete cosines
& SAB_ELASTIC_EXACT = 4 ! Exact treatment for coherent elastic
SAB_ELASTIC_DISCRETE = 3, & ! Sample from discrete cosines
SAB_ELASTIC_EXACT = 4 ! Exact treatment for coherent elastic
! Reaction types
integer, parameter :: &
@ -201,6 +208,15 @@ module constants
ASCII = 1, & ! ASCII cross section file
BINARY = 2 ! Binary cross section file
! Probability table parameters
integer, parameter :: &
URR_CUM_PROB = 1, &
URR_TOTAL = 2, &
URR_ELASTIC = 3, &
URR_FISSION = 4, &
URR_N_GAMMA = 5, &
URR_HEATING = 6
! Maximum number of partial fission reactions
integer, parameter :: PARTIAL_FISSION_MAX = 4
@ -280,10 +296,10 @@ module constants
PROB_CRITICALITY = 2 ! Criticality problem
! Unit numbers
integer, parameter :: UNIT_LOG = 11 ! unit # for writing log file
integer, parameter :: UNIT_TALLY = 12 ! unit # for writing tally file
integer, parameter :: UNIT_PLOT = 13 ! unit # for writing plot file
integer, parameter :: UNIT_CMFD = 14 ! unit # for writing cmfd file
integer, parameter :: UNIT_SUMMARY = 11 ! unit # for writing summary file
integer, parameter :: UNIT_TALLY = 12 ! unit # for writing tally file
integer, parameter :: UNIT_PLOT = 13 ! unit # for writing plot file
integer, parameter :: UNIT_CMFD = 14 ! unit # for writing cmfd file
end module constants

1572
src/cross_section.F90
View file

File diff suppressed because it is too large Load diff

16
src/endf.F90
View file

@ -1,7 +1,7 @@
module endf
use constants
use string, only: int_to_str
use string, only: to_str
contains
@ -72,7 +72,7 @@ contains
case (N_NPA)
string = '(n,npa)'
case (N_N1 : N_N40)
string = '(n,n' // trim(int_to_str(MT-50)) // ')'
string = '(n,n' // trim(to_str(MT-50)) // ')'
case (N_NC)
string = '(n,nc)'
case (N_GAMMA)
@ -122,27 +122,27 @@ contains
case (444)
string = '(damage)'
case (600 : 648)
string = '(n,p' // trim(int_to_str(MT-600)) // ')'
string = '(n,p' // trim(to_str(MT-600)) // ')'
case (649)
string = '(n,pc)'
case (650 : 698)
string = '(n,d' // trim(int_to_str(MT-650)) // ')'
string = '(n,d' // trim(to_str(MT-650)) // ')'
case (699)
string = '(n,dc)'
case (700 : 748)
string = '(n,t' // trim(int_to_str(MT-700)) // ')'
string = '(n,t' // trim(to_str(MT-700)) // ')'
case (749)
string = '(n,tc)'
case (750 : 798)
string = '(n,3He' // trim(int_to_str(MT-750)) // ')'
string = '(n,3He' // trim(to_str(MT-750)) // ')'
case (799)
string = '(n,3Hec)'
case (800 : 848)
string = '(n,a' // trim(int_to_str(MT-800)) // ')'
string = '(n,a' // trim(to_str(MT-800)) // ')'
case (849)
string = '(n,tc)'
case default
string = 'MT=' // trim(int_to_str(MT))
string = 'MT=' // trim(to_str(MT))
end select
end function reaction_name

2
src/energy_grid.F90
View file

@ -99,7 +99,7 @@ contains
! if the original list is empty, we need to allocate the first element and
! store first energy point
if (list_size(list) == 0) then
if (.not. associated(list)) then
allocate(list)
current => list
do i = 1, n

53
src/error.F90
View file

@ -2,7 +2,11 @@ module error
use ISO_FORTRAN_ENV
use global, only: master, free_memory, message
use global, only: master, free_memory, message, mpi_err
#ifdef MPI
use mpi
#endif
implicit none
@ -25,14 +29,14 @@ contains
! Only allow master to print to screen
if (.not. master) return
write(ou, fmt='(1X,A9)', advance='no') 'WARNING: '
write(OUTPUT_UNIT, fmt='(1X,A9)', advance='no') 'WARNING: '
n_lines = (len_trim(message)-1)/70 + 1
do i = 1, n_lines
if (i == 1) then
write(ou, fmt='(A70)') message(70*(i-1)+1:70*i)
write(OUTPUT_UNIT, fmt='(A70)') message(70*(i-1)+1:70*i)
else
write(ou, fmt='(10X,A70)') message(70*(i-1)+1:70*i)
write(OUTPUT_UNIT, fmt='(10X,A70)') message(70*(i-1)+1:70*i)
end if
end do
@ -44,29 +48,42 @@ contains
! the program is aborted.
!===============================================================================
subroutine fatal_error()
subroutine fatal_error(error_code)
integer, optional :: error_code ! error code
integer :: code ! error code
integer :: n_lines ! number of lines
integer :: i ! loop index over lines
! Only allow master to print to screen
if (master) then
write(eu, fmt='(1X,A7)', advance='no') 'ERROR: '
n_lines = (len_trim(message)-1)/72 + 1
do i = 1, n_lines
if (i == 1) then
write(eu, fmt='(A72)') message(72*(i-1)+1:72*i)
else
write(eu, fmt='(7X,A72)') message(72*(i-1)+1:72*i)
end if
end do
write(eu,*)
! set default error code
if (present(error_code)) then
code = error_code
else
code = -1
end if
! Only allow master to print to screen
write(ERROR_UNIT, fmt='(1X,A7)', advance='no') 'ERROR: '
n_lines = (len_trim(message)-1)/72 + 1
do i = 1, n_lines
if (i == 1) then
write(ERROR_UNIT, fmt='(A72)') message(72*(i-1)+1:72*i)
else
write(ERROR_UNIT, fmt='(7X,A72)') message(72*(i-1)+1:72*i)
end if
end do
write(ERROR_UNIT,*)
! Release memory from all allocatable arrays
call free_memory()
#ifdef MPI
! Abort MPI
call MPI_ABORT(MPI_COMM_WORLD, code, mpi_err)
#endif
! Abort program
stop

90
src/fileio.F90
View file

@ -1,90 +0,0 @@
module fileio
use constants
use global, only: MAX_LINE_LEN, MAX_WORD_LEN, MAX_WORDS
use string, only: split_string_wl, lower_case
implicit none
contains
!===============================================================================
! READ_LINE reads a line from a file open on a unit
!===============================================================================
subroutine read_line(unit, line, ioError)
integer, intent(in) :: unit ! unit to read from
character(*), intent(out) :: line ! line to return
integer, intent(out) :: ioError ! error status
read(UNIT=unit, FMT='(A)', IOSTAT=ioError) line
end subroutine read_line
!===============================================================================
! GET_NEXT_LINE reads the next line to the file connected on the specified unit
! including any continuation lines. If a line ends in an ampersand, the next
! line is read and its words are appended to the final array
!===============================================================================
subroutine get_next_line(unit, words, n, ioError)
integer, intent(in) :: unit ! unit to read from
character(*), intent(out) :: words(MAX_WORDS) ! words read
integer, intent(out) :: n ! number of words
integer, intent(out) :: ioError ! error status
character(MAX_LINE_LEN) :: line ! single line
character(MAX_WORD_LEN) :: local_words(MAX_WORDS) ! words on one line
integer :: index_word ! index of words
index_word = 0
do
! read line from file
read(UNIT=unit, FMT='(A100)', IOSTAT=ioError) line
! if we're at the end of the file, return
if (ioError /= 0) return
! split a single line into words
call split_string_wl(line, local_words, n)
! if there are no words, we're done
if (n == 0) exit
! Check whether there is a continuation line
if (local_words(n) == '&') then
words(index_word+1:index_word+n-1) = local_words(1:n-1)
index_word = index_word + n - 1
else
words(index_word+1:index_word+n) = local_words(1:n)
index_word = index_word + n
exit
end if
end do
! set total number of words
n = index_word
end subroutine get_next_line
!===============================================================================
! SKIP_LINES skips 'n_lines' lines from a file open on a unit
!===============================================================================
subroutine skip_lines(unit, n_lines, ioError)
integer, intent(in) :: unit ! unit to read from
integer, intent(in) :: n_lines ! number of lines to skip
integer, intent(out) :: ioError ! error status
integer :: i ! index for number of lines
do i = 1, n_lines
read(UNIT=unit, FMT=*, IOSTAT=ioError)
end do
end subroutine skip_lines
end module fileio

10
src/fission.F90
View file

@ -1,11 +1,11 @@
module fission
use ace_header, only: Nuclide
use constants
use cross_section_header, only: Nuclide
use error, only: fatal_error
use global, only: message
use interpolation, only: interpolate_tab1
use search, only: binary_search
use error, only: fatal_error
use global, only: message
use interpolation, only: interpolate_tab1
use search, only: binary_search
implicit none

210
src/geometry.F90
View file

@ -7,7 +7,7 @@ module geometry
use global
use output, only: write_message
use particle_header, only: Particle, LocalCoord, deallocate_coord
use string, only: int_to_str
use string, only: to_str
use tally, only: score_surface_current
implicit none
@ -171,26 +171,42 @@ contains
x = ceiling((xyz(1) - lat % x0)/lat % width_x)
y = ceiling((xyz(2) - lat % y0)/lat % width_y)
! Create new level of coordinates
p % in_lower_universe = .true.
allocate(p % coord % next)
! adjust local position of particle
p % coord % next % xyz(1) = p % coord % xyz(1) - &
(lat%x0 + (x-0.5_8)*lat%width_x)
p % coord % next % xyz(2) = p % coord % xyz(2) - &
(lat%y0 + (y-0.5_8)*lat%width_y)
p % coord % next % xyz(3) = p % coord % xyz(3)
p % coord % next % uvw = p % coord % uvw
! Check if lattice coordinates are within bounds
if (x < 1 .or. x > lat % n_x .or. &
y < 1 .or. y > lat % n_y) then
! Move particle to next level
p % coord => p % coord % next
! This condition should only get hit in rare circumstances where
! a neutron hits the corner of a lattice. In this case, the
! neutron may need to be moved diagonally across the lattice. To
! do so, we remove all lower coordinate levels and then search
! from universe 0.
p % coord => p % coord0
call deallocate_coord(p % coord % next)
else
! Create new level of coordinates
p % in_lower_universe = .true.
allocate(p % coord % next)
! set particle lattice indices
p % coord % lattice = c % fill
p % coord % lattice_x = x
p % coord % lattice_y = y
p % coord % universe = lat % element(x,y)
! adjust local position of particle
p % coord % next % xyz(1) = p % coord % xyz(1) - &
(lat%x0 + (x-0.5_8)*lat%width_x)
p % coord % next % xyz(2) = p % coord % xyz(2) - &
(lat%y0 + (y-0.5_8)*lat%width_y)
p % coord % next % xyz(3) = p % coord % xyz(3)
p % coord % next % uvw = p % coord % uvw
! Move particle to next level
p % coord => p % coord % next
! set particle lattice indices
p % coord % lattice = c % fill
p % coord % lattice_x = x
p % coord % lattice_y = y
p % coord % universe = lat % element(x,y)
end if
call find_cell(p, found)
if (.not. found) exit
@ -216,24 +232,24 @@ contains
type(Particle), pointer :: p
integer, intent(in) :: last_cell ! last cell particle was in
real(8) :: x ! x-x0 for sphere
real(8) :: y ! y-y0 for sphere
real(8) :: z ! z-z0 for sphere
real(8) :: R ! radius of sphere
real(8) :: u ! x-component of direction
real(8) :: v ! y-component of direction
real(8) :: w ! z-component of direction
real(8) :: n1 ! x-component of surface normal
real(8) :: n2 ! y-component of surface normal
real(8) :: n3 ! z-component of surface normal
real(8) :: dot_prod ! dot product of direction and normal
real(8) :: norm ! "norm" of surface normal
logical :: found ! particle found in universe?
real(8) :: x ! x-x0 for sphere
real(8) :: y ! y-y0 for sphere
real(8) :: z ! z-z0 for sphere
real(8) :: R ! radius of sphere
real(8) :: u ! x-component of direction
real(8) :: v ! y-component of direction
real(8) :: w ! z-component of direction
real(8) :: n1 ! x-component of surface normal
real(8) :: n2 ! y-component of surface normal
real(8) :: n3 ! z-component of surface normal
real(8) :: dot_prod ! dot product of direction and normal
real(8) :: norm ! "norm" of surface normal
logical :: found ! particle found in universe?
type(Surface), pointer :: surf => null()
surf => surfaces(abs(p % surface))
if (verbosity >= 10 .or. trace) then
message = " Crossing surface " // trim(int_to_str(surf % id))
message = " Crossing surface " // trim(to_str(surf % id))
call write_message()
end if
@ -258,7 +274,7 @@ contains
! Display message
if (verbosity >= 10 .or. trace) then
message = " Leaked out of surface " // trim(int_to_str(surf % id))
message = " Leaked out of surface " // trim(to_str(surf % id))
call write_message()
end if
return
@ -269,7 +285,8 @@ contains
! Do not handle reflective boundary conditions on lower universes
if (p % in_lower_universe) then
message = "Cannot reflect particle off surface in a lower universe."
message = "Cannot reflect particle " // trim(to_str(p % id)) // &
" off surface in a lower universe."
call fatal_error()
end if
@ -367,7 +384,7 @@ contains
p % coord0 % uvw = (/ u, v, w /)
case default
message = "Reflection not supported for surface " // &
trim(int_to_str(surf % id))
trim(to_str(surf % id))
call fatal_error()
end select
@ -380,7 +397,7 @@ contains
! Diagnostic message
if (verbosity >= 10 .or. trace) then
message = " Reflected from surface " // trim(int_to_str(surf%id))
message = " Reflected from surface " // trim(to_str(surf%id))
call write_message()
end if
return
@ -412,9 +429,9 @@ contains
! Couldn't find next cell anywhere!
if ((.not. found) .and. (.not. plotting)) then
message = "After particle crossed surface " // trim(int_to_str( &
surfaces(abs(p%surface)) % id)) // " it could not be located in " &
// "any cell and it did not leak."
message = "After particle " // trim(to_str(p % id)) // " crossed surface " &
// trim(to_str(surfaces(abs(p%surface)) % id)) // " it could not be &
&located in any cell and it did not leak."
call fatal_error()
end if
@ -424,110 +441,69 @@ contains
! CROSS_LATTICE moves a particle into a new lattice element
!===============================================================================
subroutine cross_lattice(p)
subroutine cross_lattice(p, lattice_crossed)
type(Particle), pointer :: p
integer, intent(in) :: lattice_crossed
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?
integer :: i_x ! x index in lattice
integer :: i_y ! y index in lattice
real(8) :: x0 ! half the width of lattice element
real(8) :: y0 ! half the height of lattice element
logical :: found ! particle found in cell?
type(Lattice), pointer :: lat => null()
lat => lattices(p % coord % lattice)
if (verbosity >= 10 .or. trace) then
message = " Crossing lattice " // trim(int_to_str(lat % id)) // &
". Current position (" // trim(int_to_str(p % coord % lattice_x)) &
// "," // trim(int_to_str(p % coord % lattice_y)) // ")"
message = " Crossing lattice " // trim(to_str(lat % id)) // &
". Current position (" // trim(to_str(p % coord % lattice_x)) &
// "," // trim(to_str(p % coord % lattice_y)) // ")"
call write_message()
end if
u = p % coord % uvw(1)
v = p % coord % uvw(2)
if (lat % type == LATTICE_RECT) then
x = p % coord % xyz(1)
y = p % coord % xyz(2)
z = p % coord % xyz(3)
x0 = lat % width_x * 0.5_8
y0 = lat % width_y * 0.5_8
dist = INFINITY
! left and right sides
if (u == ZERO) then
d_left = INFINITY
d_right = INFINITY
elseif (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 == ZERO) then
d_bottom = INFINITY
d_top = INFINITY
elseif (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
select case (lattice_crossed)
case (LATTICE_LEFT)
! Move particle to left element
p % coord % lattice_x = p % coord % lattice_x - 1
p % coord % xyz(1) = x0
elseif (dist == d_right) then
case (LATTICE_RIGHT)
! Move particle to right element
p % coord % lattice_x = p % coord % lattice_x + 1
p % coord % xyz(1) = -x0
elseif (dist == d_bottom) then
case (LATTICE_BOTTOM)
! Move particle to bottom element
p % coord % lattice_y = p % coord % lattice_y - 1
p % coord % xyz(2) = y0
elseif (dist == d_top) then
case (LATTICE_TOP)
! Move particle to top element
p % coord % lattice_y = p % coord % lattice_y + 1
p % coord % xyz(2) = -y0
end if
end select
elseif (lat % type == LATTICE_HEX) then
! TODO: Add hex lattice support
end if
! Check to make sure still in lattice
i_x = p % coord % lattice_x
i_y = p % coord % lattice_y
if (i_x < 1 .or. i_x > lat % n_x) then
message = "Reached edge of lattice " // trim(int_to_str(lat % id)) // &
" at position (" // trim(int_to_str(i_x)) // "," // &
trim(int_to_str(i_y)) // ")."
message = "Particle " // trim(to_str(p % id)) // " reached edge of &
&lattice " // trim(to_str(lat % id)) // " at position (" // &
trim(to_str(i_x)) // "," // trim(to_str(i_y)) // ")."
call fatal_error()
elseif (i_y < 1 .or. i_y > lat % n_y) then
message = "Reached edge of lattice " // trim(int_to_str(lat % id)) // &
" at position (" // trim(int_to_str(i_x)) // "," // &
trim(int_to_str(i_y)) // ")."
message = "Particle " // trim(to_str(p % id)) // " reached edge of &
&lattice " // trim(to_str(lat % id)) // " at position (" // &
trim(to_str(i_x)) // "," // trim(to_str(i_y)) // ")."
call fatal_error()
end if
@ -537,8 +513,8 @@ contains
! Find cell in next lattice element
call find_cell(p, found)
if (.not. found) then
message = "Could not locate particle in universe: " // &
int_to_str(universes(p % coord % universe) % id)
message = "Could not locate particle " // trim(to_str(p % id)) // &
" in universe " // to_str(universes(p % coord % universe) % id)
call fatal_error()
end if
@ -555,7 +531,7 @@ contains
type(Particle), pointer :: p
real(8), intent(out) :: dist
integer, intent(out) :: surface_crossed
logical, intent(out) :: lattice_crossed
integer, intent(out) :: lattice_crossed
integer :: i ! index for surface in cell
integer :: index_surf ! index in surfaces array (with sign)
@ -576,7 +552,7 @@ contains
! inialize distance to infinity (huge)
dist = INFINITY
lattice_crossed = .false.
lattice_crossed = NONE
nullify(final_coord)
! Get pointer to top-level coordinates
@ -866,7 +842,7 @@ contains
if (d < dist) then
dist = d
surface_crossed = -cl % surfaces(i)
lattice_crossed = .false.
lattice_crossed = NONE
final_coord => coord
end if
@ -906,7 +882,11 @@ contains
if (d < dist) then
if (abs(d - dist)/dist >= FP_PRECISION) then
dist = d
lattice_crossed = .true.
if (u > 0) then
lattice_crossed = LATTICE_RIGHT
else
lattice_crossed = LATTICE_LEFT
end if
final_coord => coord
end if
end if
@ -923,7 +903,11 @@ contains
if (d < dist) then
if (abs(d - dist)/dist >= FP_PRECISION) then
dist = d
lattice_crossed = .true.
if (v > 0) then
lattice_crossed = LATTICE_TOP
else
lattice_crossed = LATTICE_BOTTOM
end if
final_coord => coord
end if
end if
@ -1072,7 +1056,7 @@ contains
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
z >= z0 .and. z < z1) then
s = SENSE_NEGATIVE
else
s = SENSE_POSITIVE
@ -1091,7 +1075,7 @@ contains
I = surf % coeffs(9)
J = surf % coeffs(10)
func = A*x*x + B*y*y + C*z*z + D*x*y + E*y*z + F*x*z + G*x &
& + H*y + I*z + J
+ H*y + I*z + J
end select

4
src/geometry_header.F90
View file

@ -45,8 +45,8 @@ module geometry_header
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
neighbor_pos(:), & ! List of cells on positive side
neighbor_neg(:) ! List of cells on negative side
integer :: bc ! Boundary condition
end type Surface

38
src/global.F90
View file

@ -1,18 +1,18 @@
module global
use bank_header, only: Bank
use ace_header, only: Nuclide, SAB_Table, xsListing, NuclideMicroXS, &
MaterialMacroXS
use bank_header, only: Bank
use cmfd_header
use constants
use cross_section_header, only: Nuclide, SAB_Table, xsListing, &
NuclideMicroXS, MaterialMacroXS
use datatypes_header, only: DictionaryII, DictionaryCI
use geometry_header, only: Cell, Universe, Lattice, Surface
use material_header, only: Material
use mesh_header, only: StructuredMesh
use particle_header, only: Particle
use source_header, only: ExtSource
use tally_header, only: TallyObject, TallyMap
use timing, only: Timer
use datatypes_header, only: DictionaryII, DictionaryCI
use geometry_header, only: Cell, Universe, Lattice, Surface
use material_header, only: Material
use mesh_header, only: StructuredMesh
use particle_header, only: Particle
use source_header, only: ExtSource
use tally_header, only: TallyObject, TallyMap
use timing, only: Timer
#ifdef MPI
use mpi
@ -74,6 +74,9 @@ module global
integer :: n_grid ! number of points on unionized grid
real(8), allocatable :: e_grid(:) ! energies on unionized grid
! Unreoslved resonance probablity tables
logical :: urr_ptables_on = .false.
! ============================================================================
! TALLY-RELATED VARIABLES
@ -112,6 +115,13 @@ module global
real(8) :: keff = ONE
real(8) :: keff_std
! Shannon entropy
logical :: entropy_on = .false.
real(8) :: entropy ! value of shannon entropy
real(8) :: entropy_lower_left(3) ! lower-left corner for entropy box
real(8) :: entropy_upper_right(3) ! upper-right corner for entropy box
real(8), allocatable :: entropy_p(:,:,:)
! ============================================================================
! PARALLEL PROCESSING VARIABLES
@ -120,6 +130,7 @@ module global
logical :: master ! master process?
logical :: mpi_enabled ! is MPI in use and initialized?
integer :: mpi_err ! MPI error code
integer :: MPI_BANK ! MPI datatype for fission bank
! ============================================================================
! TIMING VARIABLES
@ -212,11 +223,6 @@ contains
if (allocated(fission_bank)) deallocate(fission_bank)
if (allocated(source_bank)) deallocate(source_bank)
#ifdef MPI
! If MPI is in use and enabled, terminate it
call MPI_FINALIZE(mpi_err)
#endif
end subroutine free_memory
end module global

142
src/initialize.F90
View file

@ -1,7 +1,8 @@
module initialize
use ace, only: read_xs
use bank_header, only: Bank
use constants
use cross_section, only: read_xs
use datatypes, only: dict_create, dict_add_key, dict_get_key, &
dict_has_key, dict_keys
use datatypes_header, only: ListKeyValueII, DictionaryII
@ -12,16 +13,18 @@ module initialize
use global
use input_xml, only: read_input_xml, read_cross_sections_xml, &
cells_in_univ_dict
use logging, only: create_log
use mpi_routines, only: setup_mpi
use output, only: title, header, print_summary, print_geometry, &
print_plot
print_plot, create_summary_file
use random_lcg, only: initialize_prng
use source, only: initialize_source
use string, only: int_to_str, starts_with, ends_with, lower_case
use string, only: to_str, starts_with, ends_with, lower_case
use tally, only: create_tally_map, TallyObject
use timing, only: timer_start, timer_stop
#ifdef MPI
use mpi
#endif
implicit none
type(DictionaryII), pointer :: build_dict => null()
@ -47,13 +50,13 @@ contains
! Read command line arguments
call read_command_line()
if (master) call create_log()
if (master) call create_summary_file()
! Print the OpenMC title and version/date/time information
if (master) call title()
! Print initialization header block
if (master) call header("INITIALIZATION", 1)
if (master) call header("INITIALIZATION", level=1)
! Initialize random number generator
call initialize_prng()
@ -120,6 +123,81 @@ contains
end subroutine initialize_run
!===============================================================================
! SETUP_MPI initilizes the Message Passing Interface (MPI) and determines the
! number of processors the problem is being run with as well as the rank of each
! processor.
!===============================================================================
subroutine setup_mpi()
#ifdef MPI
integer :: i
integer :: bank_blocks(4) ! Count for each datatype
integer :: bank_types(4) ! Datatypes
integer(MPI_ADDRESS_KIND) :: bank_disp(4) ! Displacements
integer(MPI_ADDRESS_KIND) :: base
type(Bank) :: b
mpi_enabled = .true.
! Initialize MPI
call MPI_INIT(mpi_err)
if (mpi_err /= MPI_SUCCESS) then
message = "Failed to initialize MPI."
call fatal_error()
end if
! Determine number of processors
call MPI_COMM_SIZE(MPI_COMM_WORLD, n_procs, mpi_err)
if (mpi_err /= MPI_SUCCESS) then
message = "Could not determine number of processors."
call fatal_error()
end if
! Determine rank of each processor
call MPI_COMM_RANK(MPI_COMM_WORLD, rank, mpi_err)
if (mpi_err /= MPI_SUCCESS) then
message = "Could not determine MPI rank."
call fatal_error()
end if
! Determine master
if (rank == 0) then
master = .true.
else
master = .false.
end if
! Determine displacements for MPI_BANK type
call MPI_GET_ADDRESS(b % id, bank_disp(1), mpi_err)
call MPI_GET_ADDRESS(b % xyz, bank_disp(2), mpi_err)
call MPI_GET_ADDRESS(b % uvw, bank_disp(3), mpi_err)
call MPI_GET_ADDRESS(b % E, bank_disp(4), mpi_err)
! Adjust displacements
base = bank_disp(1)
do i = 1, 4
bank_disp(i) = bank_disp(i) - base
end do
! Define MPI_BANK for fission sites
bank_blocks = (/ 1, 3, 3, 1 /)
bank_types = (/ MPI_INTEGER8, MPI_REAL8, MPI_REAL8, MPI_REAL8 /)
call MPI_TYPE_CREATE_STRUCT(4, bank_blocks, bank_disp, &
bank_types, MPI_BANK, mpi_err)
call MPI_TYPE_COMMIT(MPI_BANK, mpi_err)
#else
! if no MPI, set processor to master
mpi_enabled = .false.
rank = 0
n_procs = 1
master = .true.
#endif
end subroutine setup_mpi
!===============================================================================
! READ_COMMAND_LINE reads all parameters from the command line
!===============================================================================
@ -297,8 +375,8 @@ contains
i_array = dict_get_key(surface_dict, abs(id))
c % surfaces(j) = sign(i_array, id)
else
message = "Could not find surface " // trim(int_to_str(abs(id))) // &
& " specified on cell " // trim(int_to_str(c % id))
message = "Could not find surface " // trim(to_str(abs(id))) // &
" specified on cell " // trim(to_str(c % id))
call fatal_error()
end if
end if
@ -311,8 +389,8 @@ contains
if (dict_has_key(universe_dict, id)) then
c % universe = dict_get_key(universe_dict, id)
else
message = "Could not find universe " // trim(int_to_str(id)) // &
" specified on cell " // trim(int_to_str(c % id))
message = "Could not find universe " // trim(to_str(id)) // &
" specified on cell " // trim(to_str(c % id))
call fatal_error()
end if
@ -325,8 +403,8 @@ contains
c % type = CELL_NORMAL
c % material = dict_get_key(material_dict, id)
else
message = "Could not find material " // trim(int_to_str(id)) // &
" specified on cell " // trim(int_to_str(c % id))
message = "Could not find material " // trim(to_str(id)) // &
" specified on cell " // trim(to_str(c % id))
call fatal_error()
end if
else
@ -338,8 +416,8 @@ contains
c % type = CELL_LATTICE
c % fill = dict_get_key(lattice_dict, id)
else
message = "Specified fill " // trim(int_to_str(id)) // " on cell " // &
trim(int_to_str(c % id)) // " is neither a universe nor a lattice."
message = "Specified fill " // trim(to_str(id)) // " on cell " // &
trim(to_str(c % id)) // " is neither a universe nor a lattice."
call fatal_error()
end if
end if
@ -356,8 +434,8 @@ contains
if (dict_has_key(universe_dict, id)) then
l % element(j,k) = dict_get_key(universe_dict, id)
else
message = "Invalid universe number " // trim(int_to_str(id)) &
// " specified on lattice " // trim(int_to_str(l % id))
message = "Invalid universe number " // trim(to_str(id)) &
// " specified on lattice " // trim(to_str(l % id))
call fatal_error()
end if
end do
@ -376,8 +454,8 @@ contains
if (dict_has_key(cell_dict, id)) then
t % cell_bins(j) % scalar = dict_get_key(cell_dict, id)
else
message = "Could not find cell " // trim(int_to_str(id)) // &
& " specified on tally " // trim(int_to_str(t % id))
message = "Could not find cell " // trim(to_str(id)) // &
" specified on tally " // trim(to_str(t % id))
call fatal_error()
end if
end do
@ -392,8 +470,8 @@ contains
if (dict_has_key(surface_dict, id)) then
t % surface_bins(j) % scalar = dict_get_key(surface_dict, id)
else
message = "Could not find surface " // trim(int_to_str(id)) // &
& " specified on tally " // trim(int_to_str(t % id))
message = "Could not find surface " // trim(to_str(id)) // &
" specified on tally " // trim(to_str(t % id))
call fatal_error()
end if
end do
@ -408,8 +486,8 @@ contains
if (dict_has_key(universe_dict, id)) then
t % universe_bins(j) % scalar = dict_get_key(universe_dict, id)
else
message = "Could not find universe " // trim(int_to_str(id)) // &
& " specified on tally " // trim(int_to_str(t % id))
message = "Could not find universe " // trim(to_str(id)) // &
" specified on tally " // trim(to_str(t % id))
call fatal_error()
end if
end do
@ -424,8 +502,8 @@ contains
if (dict_has_key(material_dict, id)) then
t % material_bins(j) % scalar = dict_get_key(material_dict, id)
else
message = "Could not find material " // trim(int_to_str(id)) // &
& " specified on tally " // trim(int_to_str(t % id))
message = "Could not find material " // trim(to_str(id)) // &
" specified on tally " // trim(to_str(t % id))
call fatal_error()
end if
end do
@ -440,8 +518,8 @@ contains
if (dict_has_key(cell_dict, id)) then
t % cellborn_bins(j) % scalar = dict_get_key(cell_dict, id)
else
message = "Could not find material " // trim(int_to_str(id)) // &
& " specified on tally " // trim(int_to_str(t % id))
message = "Could not find material " // trim(to_str(id)) // &
" specified on tally " // trim(to_str(t % id))
call fatal_error()
end if
end do
@ -455,8 +533,8 @@ contains
if (dict_has_key(mesh_dict, id)) then
t % mesh = dict_get_key(mesh_dict, id)
else
message = "Could not find mesh " // trim(int_to_str(id)) // &
& " specified on tally " // trim(int_to_str(t % id))
message = "Could not find mesh " // trim(to_str(id)) // &
" specified on tally " // trim(to_str(t % id))
call fatal_error()
end if
end if
@ -546,9 +624,9 @@ contains
! Check to make sure either all atom percents or all weight percents are
! given
if (.not. (all(mat%atom_percent > ZERO) .or. &
& all(mat%atom_percent < ZERO))) then
all(mat%atom_percent < ZERO))) then
message = "Cannot mix atom and weight percents in material " // &
& int_to_str(mat % id)
to_str(mat % id)
call fatal_error()
end if
@ -606,7 +684,7 @@ contains
end do
sum_percent = ONE / sum_percent
mat % density = -mat % density * N_AVOGADRO &
& / MASS_NEUTRON * sum_percent
/ MASS_NEUTRON * sum_percent
end if
! Calculate nuclide atom densities and deallocate atom_percent array

90
src/input_xml.F90
View file

@ -9,14 +9,14 @@ module input_xml
use global
use mesh_header, only: StructuredMesh
use output, only: write_message
use string, only: lower_case, int_to_str, str_to_int, str_to_real, &
use string, only: lower_case, to_str, str_to_int, str_to_real, &
split_string, starts_with, ends_with
use tally_header, only: TallyObject
implicit none
type(DictionaryII), pointer :: & ! used to count how many cells each
& cells_in_univ_dict => null() ! universe contains
type(DictionaryII), pointer :: & ! used to count how many cells each
cells_in_univ_dict => null() ! universe contains
contains
@ -80,9 +80,8 @@ contains
! variable
call get_environment_variable("CROSS_SECTIONS", env_variable)
if (len_trim(env_variable) == 0) then
message = "No cross_sections.xml file was specified in " // &
"settings.xml or in the CROSS_SECTIONS environment " // &
"variable."
message = "No cross_sections.xml file was specified in settings.xml &
&or in the CROSS_SECTIONS environment variable."
call fatal_error()
else
path_cross_sections = trim(env_variable)
@ -144,6 +143,30 @@ contains
trace_particle = trace_(2)
end if
! Entropy box
if (associated(entropy_box_)) then
! Check to make sure enough values were supplied
if (size(entropy_box_) /= 6) then
message = "Need to supply lower-left and upper-right coordinates &
&for Shannon entropy box."
call fatal_error()
end if
! Copy values
entropy_lower_left = entropy_box_(1:3)
entropy_upper_right = entropy_box_(4:6)
! Check on values provided
if (.not. all(entropy_upper_right > entropy_lower_left)) then
message = "Upper-right coordinate must be greater than lower-left &
&coordinate for Shannon entropy box."
call fatal_error()
end if
! Turn on Shannon entropy calculation
entropy_on = .true.
end if
end subroutine read_settings_xml
!===============================================================================
@ -203,7 +226,7 @@ contains
! Check to make sure that either material or fill was specified
if (c % material == 0 .and. c % fill == 0) then
message = "Neither material nor fill was specified for cell " // &
trim(int_to_str(c % id))
trim(to_str(c % id))
call fatal_error()
end if
@ -217,7 +240,7 @@ contains
! Check to make sure that surfaces were specified
if (.not. associated(cell_(i) % surfaces)) then
message = "No surfaces specified for cell " // &
trim(int_to_str(c % id))
trim(to_str(c % id))
call fatal_error()
end if
@ -313,12 +336,11 @@ contains
n = size(surface_(i) % coeffs)
if (n < coeffs_reqd) then
message = "Not enough coefficients specified for surface: " // &
trim(int_to_str(s % id))
print *, n, coeffs_reqd
trim(to_str(s % id))
call fatal_error()
elseif (n > coeffs_reqd) then
message = "Too many coefficients specified for surface: " // &
trim(int_to_str(s % id))
trim(to_str(s % id))
call fatal_error()
else
allocate(s % coeffs(n))
@ -339,7 +361,7 @@ contains
s % bc = BC_PERIODIC
case default
message = "Unknown boundary condition '" // trim(word) // &
"' specified on surface " // trim(int_to_str(s % id))
"' specified on surface " // trim(to_str(s % id))
call fatal_error()
end select
@ -487,14 +509,14 @@ contains
m % density = 1.0e-24 * val
case default
message = "Unkwown units '" // trim(material_(i) % density % units) &
// "' specified on material " // trim(int_to_str(m % id))
// "' specified on material " // trim(to_str(m % id))
call fatal_error()
end select
! Check to ensure material has at least one nuclide
if (.not. associated(material_(i) % nuclides)) then
message = "No nuclides specified on material " // &
trim(int_to_str(m % id))
trim(to_str(m % id))
call fatal_error()
end if
@ -514,7 +536,7 @@ contains
! Check for empty name on nuclide
if (len_trim(nuc % name) == 0) then
message = "No name specified on nuclide in material " // &
trim(int_to_str(m % id))
trim(to_str(m % id))
call fatal_error()
end if
@ -522,7 +544,7 @@ contains
if (len_trim(nuc % xs) == 0) then
if (default_xs == '') then
message = "No cross section specified for nuclide in material " &
// trim(int_to_str(m % id))
// trim(to_str(m % id))
call fatal_error()
else
nuc % xs = default_xs
@ -540,8 +562,8 @@ contains
trim(name)
call fatal_error()
elseif (nuc % ao /= ZERO .and. nuc % wo /= ZERO) then
message = "Cannot specify both atom and weight percents for a " &
// "nuclide: " // trim(name)
message = "Cannot specify both atom and weight percents for a &
&nuclide: " // trim(name)
call fatal_error()
end if
@ -671,16 +693,16 @@ contains
! Read mesh origin location
if (m % n_dimension /= size(mesh_(i) % origin)) then
message = "Number of entries on <origin> must be the same as " // &
"the number of entries on <dimension>."
message = "Number of entries on <origin> must be the same as the &
&number of entries on <dimension>."
call fatal_error()
end if
m % origin = mesh_(i) % origin
! Read mesh widths
if (size(mesh_(i) % width) /= size(mesh_(i) % origin)) then
message = "Number of entries on <width> must be the same as " // &
"the number of entries on <origin>."
message = "Number of entries on <width> must be the same as the &
&number of entries on <origin>."
call fatal_error()
end if
m % width = mesh_(i) % width
@ -713,7 +735,7 @@ contains
if (len_trim(tally_(i) % filters % cell) > 0 .and. &
len_trim(tally_(i) % filters % surface) > 0) then
message = "Cannot specify both cell and surface filters for tally " &
// trim(int_to_str(t % id))
// trim(to_str(t % id))
call fatal_error()
end if
@ -768,8 +790,8 @@ contains
i_mesh = dict_get_key(mesh_dict, id)
m => meshes(i_mesh)
else
message = "Could not find mesh " // trim(int_to_str(id)) // &
" specified on tally " // trim(int_to_str(t % id))
message = "Could not find mesh " // trim(to_str(id)) // &
" specified on tally " // trim(to_str(t % id))
call fatal_error()
end if
@ -823,8 +845,8 @@ contains
case ('total')
t % macro_bins(j) % scalar = MACRO_TOTAL
if (t % n_bins(T_ENERGYOUT) > 0) then
message = "Cannot tally total reaction rate with an " &
// "outgoing energy filter."
message = "Cannot tally total reaction rate with an &
&outgoing energy filter."
call fatal_error()
end if
case ('scatter')
@ -848,15 +870,15 @@ contains
case ('absorption')
t % macro_bins(j) % scalar = MACRO_ABSORPTION
if (t % n_bins(T_ENERGYOUT) > 0) then
message = "Cannot tally absorption rate with an outgoing " &
// "energy filter."
message = "Cannot tally absorption rate with an outgoing &
&energy filter."
call fatal_error()
end if
case ('fission')
t % macro_bins(j) % scalar = MACRO_FISSION
if (t % n_bins(T_ENERGYOUT) > 0) then
message = "Cannot tally fission rate with an outgoing " &
// "energy filter."
message = "Cannot tally fission rate with an outgoing &
&energy filter."
call fatal_error()
end if
case ('nu-fission')
@ -868,9 +890,9 @@ contains
! Check to make sure that current is the only desired response
! for this tally
if (n_words > 1) then
message = "Cannot tally other macro reactions in the same " &
// "tally as surface currents. Separate other macro " &
// "reactions into a distinct tally."
message = "Cannot tally other macro reactions in the same &
&tally as surface currents. Separate other macro &
&reactions into a distinct tally."
call fatal_error()
end if

291
src/mpi_routines.F90 → src/intercycle.F90
View file

@ -1,7 +1,8 @@
module mpi_routines
module intercycle
use constants, only: MAX_LINE_LEN
use error, only: fatal_error
use ISO_FORTRAN_ENV
use error, only: fatal_error, warning
use global
use output, only: write_message
use particle_header, only: Particle, initialize_particle
@ -13,88 +14,8 @@ module mpi_routines
use mpi
#endif
implicit none
integer :: MPI_BANK ! MPI datatype for fission bank
integer(8) :: bank_index ! Fission bank site unique identifier
contains
!===============================================================================
! SETUP_MPI initilizes the Message Passing Interface (MPI) and determines the
! number of processors the problem is being run with as well as the rank of each
! processor.
!===============================================================================
subroutine setup_mpi()
#ifdef MPI
integer :: i
integer :: bank_blocks(4) ! Count for each datatype
integer :: bank_types(4) ! Datatypes
integer(MPI_ADDRESS_KIND) :: bank_disp(4) ! Displacements
integer(MPI_ADDRESS_KIND) :: base
type(Bank) :: b
mpi_enabled = .true.
! Initialize MPI
call MPI_INIT(mpi_err)
if (mpi_err /= MPI_SUCCESS) then
message = "Failed to initialize MPI."
call fatal_error()
end if
! Determine number of processors
call MPI_COMM_SIZE(MPI_COMM_WORLD, n_procs, mpi_err)
if (mpi_err /= MPI_SUCCESS) then
message = "Could not determine number of processors."
call fatal_error()
end if
! Determine rank of each processor
call MPI_COMM_RANK(MPI_COMM_WORLD, rank, mpi_err)
if (mpi_err /= MPI_SUCCESS) then
message = "Could not determine MPI rank."
call fatal_error()
end if
! Determine master
if (rank == 0) then
master = .true.
else
master = .false.
end if
! Determine displacements for MPI_BANK type
call MPI_GET_ADDRESS(b % id, bank_disp(1), mpi_err)
call MPI_GET_ADDRESS(b % xyz, bank_disp(2), mpi_err)
call MPI_GET_ADDRESS(b % uvw, bank_disp(3), mpi_err)
call MPI_GET_ADDRESS(b % E, bank_disp(4), mpi_err)
! Adjust displacements
base = bank_disp(1)
do i = 1, 4
bank_disp(i) = bank_disp(i) - base
end do
! Define MPI_BANK for fission sites
bank_blocks = (/ 1, 3, 3, 1 /)
bank_types = (/ MPI_INTEGER8, MPI_REAL8, MPI_REAL8, MPI_REAL8 /)
call MPI_TYPE_CREATE_STRUCT(4, bank_blocks, bank_disp, &
& bank_types, MPI_BANK, mpi_err)
call MPI_TYPE_COMMIT(MPI_BANK, mpi_err)
#else
! if no MPI, set processor to master
mpi_enabled = .false.
rank = 0
n_procs = 1
master = .true.
#endif
end subroutine setup_mpi
!===============================================================================
! SYNCHRONIZE_BANK samples source sites from the fission sites that were
! accumulated during the cycle. This routine is what allows this Monte Carlo to
@ -115,9 +36,9 @@ contains
integer(8) :: sites_needed ! # of sites to be sampled
real(8) :: p_sample ! probability of sampling a site
type(Bank), allocatable :: &
& temp_sites(:), & ! local array of extra sites on each node
& left_bank(:), & ! bank sites to send/recv to or from left node
& right_bank(:) ! bank sites to send/recv to or fram right node
temp_sites(:), & ! local array of extra sites on each node
left_bank(:), & ! bank sites to send/recv to or from left node
right_bank(:) ! bank sites to send/recv to or fram right node
#ifdef MPI
integer :: status(MPI_STATUS_SIZE) ! message status
@ -133,11 +54,11 @@ contains
! Determine starting index for fission bank and total sites in fission bank
start = 0_8
call MPI_EXSCAN(n_bank, start, 1, MPI_INTEGER8, MPI_SUM, &
& MPI_COMM_WORLD, mpi_err)
MPI_COMM_WORLD, mpi_err)
finish = start + n_bank
total = finish
call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, &
& MPI_COMM_WORLD, mpi_err)
MPI_COMM_WORLD, mpi_err)
#else
start = 0_8
@ -200,11 +121,11 @@ contains
#ifdef MPI
start = 0_8
call MPI_EXSCAN(index_local, start, 1, MPI_INTEGER8, MPI_SUM, &
& MPI_COMM_WORLD, mpi_err)
MPI_COMM_WORLD, mpi_err)
finish = start + index_local
total = finish
call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, &
& MPI_COMM_WORLD, mpi_err)
MPI_COMM_WORLD, mpi_err)
#else
start = 0_8
finish = index_local
@ -259,18 +180,18 @@ contains
if (send_to_right > 0) then
i = index_local - send_to_right + 1
call MPI_ISEND(temp_sites(i), send_to_right, MPI_BANK, rank+1, 0, &
& MPI_COMM_WORLD, request, mpi_err)
MPI_COMM_WORLD, request, mpi_err)
else if (send_to_right < 0) then
call MPI_IRECV(right_bank, -send_to_right, MPI_BANK, rank+1, 1, &
& MPI_COMM_WORLD, request_right, mpi_err)
MPI_COMM_WORLD, request_right, mpi_err)
end if
if (send_to_left < 0) then
call MPI_IRECV(left_bank, -send_to_left, MPI_BANK, rank-1, 0, &
& MPI_COMM_WORLD, request_left, mpi_err)
MPI_COMM_WORLD, request_left, mpi_err)
else if (send_to_left > 0) then
call MPI_ISEND(temp_sites(1), send_to_left, MPI_BANK, rank-1, 1, &
& MPI_COMM_WORLD, request, mpi_err)
MPI_COMM_WORLD, request, mpi_err)
end if
#endif
@ -386,7 +307,8 @@ contains
tally_temp = t % scores(:,:) % val_history
if (master) then
! Description of MPI_IN_PLANE
! The MPI_IN_PLACE specifier allows the master to copy values into a
! receive buffer without having a temporary variable
call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
@ -408,4 +330,181 @@ contains
end subroutine reduce_tallies
#endif
end module mpi_routines
!===============================================================================
! SHANNON_ENTROPY calculates the Shannon entropy of the fission source
! distribution to assess source convergence
!===============================================================================
subroutine shannon_entropy()
integer :: i ! x-index for entropy mesh
integer :: j ! y-index for entropy mesh
integer :: k ! z-index for entropy mesh
integer :: m ! index for bank sites
integer(8) :: total_bank ! total # of fission bank sites
integer, save :: n_box ! total # of boxes on mesh
integer, save :: n ! # of boxes in each dimension
real(8), save :: width(3) ! width of box in each dimension
logical :: outside_box ! were there sites outside entropy box?
! On the first pass through this subroutine, we need to determine how big
! the entropy mesh should be in each direction and then allocate a
! three-dimensional array to store the fraction of source sites in each mesh
! box
if (.not. allocated(entropy_p)) then
! determine number of boxes in each direction
n = ceiling((n_particles/20)**(1.0/3.0))
n_box = n*n*n
! determine width
width = (entropy_upper_right - entropy_lower_left)/n
! allocate p
allocate(entropy_p(n,n,n))
end if
! initialize p
entropy_p = ZERO
outside_box = .false.
! loop over fission sites and count how many are in each mesh box
FISSION_SITES: do m = 1, int(n_bank,4)
! determine indices for entropy mesh box
i = int((fission_bank(m) % xyz(1) - entropy_lower_left(1))/width(1)) + 1
j = int((fission_bank(m) % xyz(2) - entropy_lower_left(2))/width(2)) + 1
k = int((fission_bank(m) % xyz(3) - entropy_lower_left(3))/width(3)) + 1
! if outside mesh, skip particle
if (i < 1 .or. i > n .or. j < 1 .or. &
j > n .or. k < 1 .or. k > n) then
outside_box = .true.
cycle
end if
! add to appropriate mesh box
entropy_p(i,j,k) = entropy_p(i,j,k) + 1
end do FISSION_SITES
! display warning message if there were sites outside entropy box
if (outside_box) then
message = "Fission source site(s) outside of entropy box."
call warning()
end if
#ifdef MPI
! collect values from all processors
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, entropy_p, n_box, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
else
call MPI_REDUCE(entropy_p, entropy_p, n_box, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
end if
! determine total number of bank sites
call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, &
MPI_COMM_WORLD, mpi_err)
#else
total_bank = n_bank
#endif
! sum values to obtain shannon entropy
if (master) then
entropy_p = entropy_p / total_bank
entropy = -sum(entropy_p * log(entropy_p)/log(2.0), entropy_p > ZERO)
end if
end subroutine shannon_entropy
!===============================================================================
! CALCULATE_KEFF calculates the single cycle estimate of keff as well as the
! mean and standard deviation of the mean for active cycles and displays them
!===============================================================================
subroutine calculate_keff(i_cycle)
integer, intent(in) :: i_cycle ! index of current cycle
integer(8) :: total_bank ! total number of source sites
integer :: n ! active cycle number
real(8) :: k_cycle ! single cycle estimate of keff
real(8), save :: k_sum ! accumulated keff
real(8), save :: k_sum_sq ! accumulated keff**2
message = "Calculate cycle keff..."
call write_message(8)
! initialize sum and square of sum at beginning of run
if (i_cycle == 1) then
k_sum = ZERO
k_sum_sq = ZERO
end if
#ifdef MPI
! Collect number bank sites onto master process
call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, &
MPI_COMM_WORLD, mpi_err)
#else
total_bank = n_bank
#endif
! Collect statistics and print output
if (master) then
! Since the creation of bank sites was originally weighted by the last
! cycle keff, we need to multiply by that keff to get the current cycle's
! value
k_cycle = real(total_bank)/real(n_particles)*keff
if (i_cycle > n_inactive) then
! Active cycle number
n = i_cycle - n_inactive
! Accumulate cycle estimate of k
k_sum = k_sum + k_cycle
k_sum_sq = k_sum_sq + k_cycle*k_cycle
! Determine mean and standard deviation of mean
keff = k_sum/n
keff_std = sqrt((k_sum_sq/n - keff*keff)/n)
! Display output for this cycle
if (i_cycle > n_inactive+1) then
if (entropy_on) then
write(UNIT=OUTPUT_UNIT, FMT=103) i_cycle, k_cycle, entropy, &
keff, keff_std
else
write(UNIT=OUTPUT_UNIT, FMT=101) i_cycle, k_cycle, keff, keff_std
end if
else
if (entropy_on) then
write(UNIT=OUTPUT_UNIT, FMT=102) i_cycle, k_cycle, entropy
else
write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle
end if
end if
else
! Display output for inactive cycle
if (entropy_on) then
write(UNIT=OUTPUT_UNIT, FMT=102) i_cycle, k_cycle, entropy
else
write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle
end if
keff = k_cycle
end if
end if
#ifdef MPI
! Broadcast new keff value to all processors
call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err)
#endif
100 format (2X,I5,2X,F8.5)
101 format (2X,I5,2X,F8.5,5X,F8.5," +/-",F8.5)
102 format (2X,I5,2X,F8.5,3X,F8.5)
103 format (2X,I5,2X,F8.5,3X,F8.5,3X,F8.5," +/-",F8.5)
end subroutine calculate_keff
end module intercycle

6
src/interpolation.F90
View file

@ -5,7 +5,7 @@ module interpolation
use error, only: fatal_error
use global, only: message
use search, only: binary_search
use string, only: int_to_str
use string, only: to_str
implicit none
@ -118,7 +118,7 @@ contains
r = (log(x) - log(x0))/(log(x1) - log(x0))
y = exp((1-r)*log(y0) + r*log(y1))
case default
message = "Unsupported interpolation scheme: " // int_to_str(interp)
message = "Unsupported interpolation scheme: " // to_str(interp)
call fatal_error()
end select
@ -204,7 +204,7 @@ contains
r = (log(x) - log(x0))/(log(x1) - log(x0))
y = exp((1-r)*log(y0) + r*log(y1))
case default
message = "Unsupported interpolation scheme: " // int_to_str(interp)
message = "Unsupported interpolation scheme: " // to_str(interp)
call fatal_error()
end select

43
src/logging.F90
View file

@ -1,43 +0,0 @@
module logging
use constants, only: MAX_WORD_LEN, UNIT_LOG
use global, only: path_input
implicit none
contains
!===============================================================================
! CREATE_LOG creates a new log file (or overwrites the existing log)
!===============================================================================
subroutine create_log()
character(MAX_WORD_LEN) :: path_log ! path of log file
logical :: file_exists ! does log file already exist?
! integer :: ioError ! error status for file access
! Create filename for log file
path_log = trim(path_input) // ".log"
! Check if log file already exists
inquire(FILE=path_log, EXIST=file_exists)
if (file_exists) then
! Possibly copy old log file
end if
! Open log file for writing
! open(FILE=path_log, UNIT=UNIT_LOG, STATUS='replace', &
! & ACTION='write', IOSTAT=ioError)
end subroutine create_log
!===============================================================================
! LOG_TALLIES
!===============================================================================
subroutine log_tallies()
end subroutine log_tallies
end module logging

35
src/main.F90
View file

@ -4,16 +4,16 @@ program main
use constants
use global
use initialize, only: initialize_run
use mpi_routines, only: synchronize_bank
use intercycle, only: shannon_entropy, calculate_keff, synchronize_bank
use output, only: write_message, header, print_runtime
use particle_header, only: Particle
use plot, only: run_plot
use physics, only: transport
use random_lcg, only: set_particle_seed
use source, only: get_source_particle
use string, only: int_to_str
use string, only: to_str
use tally, only: synchronize_tallies, write_tallies, &
tally_statistics, calculate_keff
tally_statistics
use timing, only: timer_start, timer_stop
#ifdef MPI
@ -48,6 +48,11 @@ program main
! deallocate arrays
call free_memory()
#ifdef MPI
! If MPI is in use and enabled, terminate it
call MPI_FINALIZE(mpi_err)
#endif
contains
@ -62,11 +67,24 @@ contains
integer(8) :: i_particle ! history index
type(Particle), pointer :: p => null()
if (master) call header("BEGIN SIMULATION", 1)
if (master) call header("BEGIN SIMULATION", level=1)
tallies_on = .false.
call timer_start(time_inactive)
! Display column titles
if (entropy_on) then
message = " Cycle k(cycle) Entropy Average k"
call write_message(1)
message = " ===== ======== ======= ==================="
call write_message(1)
else
message = " Cycle k(cycle) Average k"
call write_message(1)
message = " ===== ======== ==================="
call write_message(1)
end if
! ==========================================================================
! LOOP OVER CYCLES
CYCLE_LOOP: do i_cycle = 1, n_cycles
@ -74,9 +92,9 @@ contains
! Start timer for computation
call timer_start(time_compute)
message = "Simulating cycle " // trim(int_to_str(i_cycle)) // "..."
message = "Simulating cycle " // trim(to_str(i_cycle)) // "..."
call write_message(8)
! Set all tallies to zero
n_bank = 0
@ -121,6 +139,9 @@ contains
call timer_stop(time_ic_tallies)
end if
! Calculate shannon entropy
if (entropy_on) call shannon_entropy()
! Distribute fission bank across processors evenly
call synchronize_bank(i_cycle)
@ -146,7 +167,7 @@ contains
! ==========================================================================
! END OF RUN WRAPUP
if (master) call header("SIMULATION FINISHED", 1)
if (master) call header("SIMULATION FINISHED", level=1)
end subroutine run_problem

575
src/output.F90
View file

@ -2,16 +2,16 @@ module output
use ISO_FORTRAN_ENV
use ace_header, only: Nuclide, Reaction, UrrData
use constants
use cross_section_header, only: Reaction
use datatypes, only: dict_get_key
use endf, only: reaction_name
use geometry_header, only: Cell, Universe, Surface
use datatypes, only: dict_get_key
use endf, only: reaction_name
use geometry_header, only: Cell, Universe, Surface
use global
use mesh_header, only: StructuredMesh
use particle_header, only: Particle, LocalCoord
use string, only: upper_case, int_to_str, real_to_str
use tally_header, only: TallyObject
use mesh_header, only: StructuredMesh
use particle_header, only: Particle, LocalCoord
use string, only: upper_case, to_str
use tally_header, only: TallyObject
implicit none
@ -31,31 +31,46 @@ contains
character(10) :: today_date
character(8) :: today_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(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') &
' .d88888b. 888b d888 .d8888b.', &
' d88P" "Y88b 8888b d8888 d88P Y88b', &
' 888 888 88888b.d88888 888 888', &
' 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ', &
' 888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ', &
' 888 888 888 888 88888888 888 888 888 Y8P 888 888 888', &
' Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P', &
' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P"', &
'__________________888______________________________________________________', &
' 888', &
' 888'
! 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(UNIT=OUTPUT_UNIT, FMT=*) &
' Developed At: Massachusetts Institute of Technology'
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",7X,I1,".",I1,".",I1)') &
VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
! Write the date and time
call get_today(today_date, today_time)
write(ou,101) trim(today_date), trim(today_time)
101 format (6X,"Date/Time:",7X,A,1X,A)
write(ou,*)
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A,1X,A)') &
trim(today_date), trim(today_time)
! Write information to summary file
call header("OpenMC Monte Carlo Code", unit=UNIT_SUMMARY, level=1)
write(UNIT=UNIT_SUMMARY, FMT=*) &
"Copyright: 2011 Massachusetts Institute of Technology"
write(UNIT=UNIT_SUMMARY, FMT='(1X,A,7X,2(I1,"."),I1)') &
"Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
write(UNIT=UNIT_SUMMARY, FMT='(1X,"Date/Time:",5X,A,1X,A)') &
trim(today_date), trim(today_time)
! Write information on number of processors
#ifdef MPI
write(UNIT=OUTPUT_UNIT, FMT='(1X,A)') ' MPI Processes: ' // &
trim(to_str(n_procs))
write(UNIT=UNIT_SUMMARY, FMT='(1X,"MPI Processes:",1X,A)') &
trim(to_str(n_procs))
#endif
end subroutine title
@ -64,34 +79,32 @@ contains
! specified, it is assumed to be a minor header block (H3).
!===============================================================================
subroutine header(msg, level, unit_file)
subroutine header(msg, unit, level)
character(*), intent(in) :: msg
integer, optional :: level
integer, optional :: unit_file
integer, optional :: unit
integer, optional :: level
integer :: n
integer :: m
integer :: unit_
integer :: header_level
integer :: n, m
integer :: unit
character(75) :: line
character(MAX_LINE_LEN) :: line
! set default header level
if (.not. present(level)) then
header_level = 3
else
! set default level
if (present(level)) then
header_level = level
else
header_level = 3
end if
! set default unit
if (present(unit_file)) then
unit = unit_file
if (present(unit)) then
unit_ = unit
else
unit = ou
unit_ = OUTPUT_UNIT
end if
! Print first blank line
write(unit,*)
! determine how many times to repeat '=' character
n = (63 - len_trim(msg))/2
m = n
@ -104,23 +117,16 @@ contains
! print header based on level
select case (header_level)
case (1)
! determine number of spaces to put in from of header
write(unit,*) repeat('=', 75)
write(unit,*) repeat('=', n) // '> ' // trim(line) // ' <' // &
& repeat('=', m)
write(unit,*) repeat('=', 75)
write(UNIT=unit_, FMT='(/3(1X,A/))') repeat('=', 75), &
repeat('=', n) // '> ' // trim(line) // ' <' // &
repeat('=', m), repeat('=', 75)
case (2)
write(unit,*) trim(line)
write(unit,*) repeat('-', 75)
write(UNIT=unit_, FMT='(/2(1X,A/))') trim(line), repeat('-', 75)
case (3)
n = (63 - len_trim(line))/2
write(unit,*) repeat('=', n) // '> ' // trim(line) // ' <' // &
& repeat('=', m)
write(UNIT=unit_, FMT='(/1X,A/)') repeat('=', n) // '> ' // &
trim(line) // ' <' // repeat('=', m)
end select
! Print trailing blank line
write(unit, *)
end subroutine header
!===============================================================================
@ -207,13 +213,13 @@ contains
! display type of particle
select case (p % type)
case (NEUTRON)
write(ou,*) 'Neutron ' // int_to_str(p % id)
write(ou,*) 'Neutron ' // to_str(p % id)
case (PHOTON)
write(ou,*) 'Photon ' // int_to_str(p % id)
write(ou,*) 'Photon ' // to_str(p % id)
case (ELECTRON)
write(ou,*) 'Electron ' // int_to_str(p % id)
write(ou,*) 'Electron ' // to_str(p % id)
case default
write(ou,*) 'Unknown Particle ' // int_to_str(p % id)
write(ou,*) 'Unknown Particle ' // to_str(p % id)
end select
! loop through each level of universes
@ -221,26 +227,26 @@ contains
i = 0
do while(associated(coord))
! Print level
write(ou,*) ' Level ' // trim(int_to_str(i))
write(ou,*) ' Level ' // trim(to_str(i))
! Print cell for this level
if (coord % cell /= NONE) then
c => cells(coord % cell)
write(ou,*) ' Cell = ' // trim(int_to_str(c % id))
write(ou,*) ' Cell = ' // trim(to_str(c % id))
end if
! Print universe for this level
if (coord % universe /= NONE) then
u => universes(coord % universe)
write(ou,*) ' Universe = ' // trim(int_to_str(u % id))
write(ou,*) ' Universe = ' // trim(to_str(u % id))
end if
! Print information on lattice
if (coord % lattice /= NONE) then
l => lattices(coord % lattice)
write(ou,*) ' Lattice = ' // trim(int_to_str(l % id))
write(ou,*) ' Lattice position = (' // trim(int_to_str(&
p % coord % lattice_x)) // ',' // trim(int_to_str(&
write(ou,*) ' Lattice = ' // trim(to_str(l % id))
write(ou,*) ' Lattice position = (' // trim(to_str(&
p % coord % lattice_x)) // ',' // trim(to_str(&
p % coord % lattice_y)) // ')'
end if
@ -255,13 +261,13 @@ contains
! Print surface
if (p % surface /= NONE) then
s => surfaces(p % surface)
write(ou,*) ' Surface = ' // int_to_str(s % id)
write(ou,*) ' Surface = ' // to_str(s % id)
end if
write(ou,*) ' Weight = ' // real_to_str(p % wgt)
write(ou,*) ' Energy = ' // real_to_str(p % E)
write(ou,*) ' IE = ' // int_to_str(p % IE)
write(ou,*) ' Interpolation factor = ' // real_to_str(p % interp)
write(ou,*) ' Weight = ' // to_str(p % wgt)
write(ou,*) ' Energy = ' // to_str(p % E)
write(ou,*) ' IE = ' // to_str(p % IE)
write(ou,*) ' Interpolation factor = ' // to_str(p % interp)
write(ou,*)
end subroutine print_particle
@ -275,12 +281,12 @@ contains
type(Reaction), pointer :: rxn
write(ou,*) 'Reaction ' // reaction_name(rxn % MT)
write(ou,*) ' MT = ' // int_to_str(rxn % MT)
write(ou,*) ' Q-value = ' // real_to_str(rxn % Q_value)
write(ou,*) ' TY = ' // int_to_str(rxn % TY)
write(ou,*) ' Starting index = ' // int_to_str(rxn % IE)
write(ou,*) ' MT = ' // to_str(rxn % MT)
write(ou,*) ' Q-value = ' // to_str(rxn % Q_value)
write(ou,*) ' TY = ' // to_str(rxn % TY)
write(ou,*) ' Starting index = ' // to_str(rxn % IE)
if (rxn % has_energy_dist) then
write(ou,*) ' Energy: Law ' // int_to_str(rxn % edist % law)
write(ou,*) ' Energy: Law ' // to_str(rxn % edist % law)
end if
write(ou,*)
@ -290,39 +296,47 @@ contains
! PRINT_CELL displays the attributes of a cell
!===============================================================================
subroutine print_cell(c)
subroutine print_cell(c, unit)
type(Cell), pointer :: c
integer, optional :: unit
integer :: temp
integer :: i
integer :: temp
integer :: i
integer :: unit_
character(MAX_LINE_LEN) :: string
type(Universe), pointer :: u => null()
type(Lattice), pointer :: l => null()
type(Material), pointer :: m => null()
write(ou,*) 'Cell ' // int_to_str(c % id)
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
write(unit_,*) 'Cell ' // to_str(c % id)
temp = dict_get_key(cell_dict, c % id)
write(ou,*) ' Array Index = ' // int_to_str(temp)
write(unit_,*) ' Array Index = ' // to_str(temp)
u => universes(c % universe)
write(ou,*) ' Universe = ' // int_to_str(u % id)
write(unit_,*) ' Universe = ' // to_str(u % id)
select case (c % type)
case (CELL_NORMAL)
write(ou,*) ' Fill = NONE'
write(unit_,*) ' Fill = NONE'
case (CELL_FILL)
u => universes(c % fill)
write(ou,*) ' Fill = Universe ' // int_to_str(u % id)
write(unit_,*) ' Fill = Universe ' // to_str(u % id)
case (CELL_LATTICE)
l => lattices(c % fill)
write(ou,*) ' Fill = Lattice ' // int_to_str(l % id)
write(unit_,*) ' Fill = Lattice ' // to_str(l % id)
end select
if (c % material == 0) then
write(ou,*) ' Material = NONE'
write(unit_,*) ' Material = NONE'
else
m => materials(c % material)
write(ou,*) ' Material = ' // int_to_str(m % id)
write(unit_,*) ' Material = ' // to_str(m % id)
end if
write(ou,*) ' Parent Cell = ' // int_to_str(c % parent)
write(unit_,*) ' Parent Cell = ' // to_str(c % parent)
string = ""
do i = 1, c % n_surfaces
select case (c % surfaces(i))
@ -335,11 +349,11 @@ contains
case (OP_DIFFERENCE)
string = trim(string) // ' !'
case default
string = trim(string) // ' ' // int_to_str(c % surfaces(i))
string = trim(string) // ' ' // to_str(c % surfaces(i))
end select
end do
write(ou,*) ' Surface Specification:' // trim(string)
write(ou,*)
write(unit_,*) ' Surface Specification:' // trim(string)
write(unit_,*)
end subroutine print_cell
@ -347,23 +361,31 @@ contains
! PRINT_UNIVERSE displays the attributes of a universe
!===============================================================================
subroutine print_universe(univ)
subroutine print_universe(univ, unit)
type(Universe), pointer :: univ
integer, optional :: unit
integer :: i
integer :: i
integer :: unit_
character(MAX_LINE_LEN) :: string
type(Cell), pointer :: c => null()
write(ou,*) 'Universe ' // int_to_str(univ % id)
write(ou,*) ' Level = ' // int_to_str(univ % level)
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
write(unit_,*) 'Universe ' // to_str(univ % id)
write(unit_,*) ' Level = ' // to_str(univ % level)
string = ""
do i = 1, univ % n_cells
c => cells(univ % cells(i))
string = trim(string) // ' ' // int_to_str(c % id)
string = trim(string) // ' ' // to_str(c % id)
end do
write(ou,*) ' Cells =' // trim(string)
write(ou,*)
write(unit_,*) ' Cells =' // trim(string)
write(unit_,*)
end subroutine print_universe
@ -371,18 +393,27 @@ contains
! PRINT_LATTICE displays the attributes of a lattice
!===============================================================================
subroutine print_lattice(lat)
subroutine print_lattice(lat, unit)
type(Lattice), pointer :: lat
integer, optional :: unit
write(ou,*) 'Lattice ' // int_to_str(lat % id)
write(ou,*) ' n_x = ' // int_to_str(lat % n_x)
write(ou,*) ' n_y = ' // int_to_str(lat % n_y)
write(ou,*) ' x0 = ' // real_to_str(lat % x0)
write(ou,*) ' y0 = ' // real_to_str(lat % y0)
write(ou,*) ' width_x = ' // real_to_str(lat % width_x)
write(ou,*) ' width_y = ' // real_to_str(lat % width_y)
write(ou,*)
integer :: unit_
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
write(unit_,*) 'Lattice ' // to_str(lat % id)
write(unit_,*) ' n_x = ' // to_str(lat % n_x)
write(unit_,*) ' n_y = ' // to_str(lat % n_y)
write(unit_,*) ' x0 = ' // to_str(lat % x0)
write(unit_,*) ' y0 = ' // to_str(lat % y0)
write(unit_,*) ' width_x = ' // to_str(lat % width_x)
write(unit_,*) ' width_y = ' // to_str(lat % width_y)
write(unit_,*)
end subroutine print_lattice
@ -390,14 +421,22 @@ contains
! PRINT_SURFACE displays the attributes of a surface
!===============================================================================
subroutine print_surface(surf)
subroutine print_surface(surf, unit)
type(Surface), pointer :: surf
integer, optional :: unit
integer :: i
integer :: unit_
character(MAX_LINE_LEN) :: string
write(ou,*) 'Surface ' // int_to_str(surf % id)
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
write(unit_,*) 'Surface ' // to_str(surf % id)
select case (surf % type)
case (SURF_PX)
string = "X Plane"
@ -422,40 +461,40 @@ contains
case (SURF_GQ)
string = "General Quadratic"
end select
write(ou,*) ' Type = ' // trim(string)
write(unit_,*) ' Type = ' // trim(string)
string = ""
do i = 1, size(surf % coeffs)
string = trim(string) // ' ' // real_to_str(surf % coeffs(i), 4)
string = trim(string) // ' ' // to_str(surf % coeffs(i), 4)
end do
write(ou,*) ' Coefficients = ' // trim(string)
write(unit_,*) ' Coefficients = ' // 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))
string = trim(string) // ' ' // to_str(surf % neighbor_pos(i))
end do
end if
write(ou,*) ' Positive Neighbors = ' // trim(string)
write(unit_,*) ' 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))
string = trim(string) // ' ' // to_str(surf % neighbor_neg(i))
end do
end if
write(ou,*) ' Negative Neighbors =' // trim(string)
write(unit_,*) ' Negative Neighbors =' // trim(string)
select case (surf % bc)
case (BC_TRANSMIT)
write(ou,*) ' Boundary Condition = Transmission'
write(unit_,*) ' Boundary Condition = Transmission'
case (BC_VACUUM)
write(ou,*) ' Boundary Condition = Vacuum'
write(unit_,*) ' Boundary Condition = Vacuum'
case (BC_REFLECT)
write(ou,*) ' Boundary Condition = Reflective'
write(unit_,*) ' Boundary Condition = Reflective'
case (BC_PERIODIC)
write(ou,*) ' Boundary Condition = Periodic'
write(unit_,*) ' Boundary Condition = Periodic'
end select
write(ou,*)
write(unit_,*)
end subroutine print_surface
@ -463,37 +502,45 @@ contains
! PRINT_MATERIAL displays the attributes of a material
!===============================================================================
subroutine print_material(mat)
subroutine print_material(mat, unit)
type(Material), pointer :: mat
integer, optional :: unit
integer :: i
real(8) :: density
integer :: i
integer :: unit_
real(8) :: density
character(MAX_LINE_LEN) :: string
type(Nuclide), pointer :: nuc => null()
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Write identifier for material
write(ou,*) 'Material ' // int_to_str(mat % id)
write(unit_,*) 'Material ' // to_str(mat % id)
! Write total atom density in atom/b-cm
write(ou,*) ' Atom Density = ' // trim(real_to_str(mat % density)) &
& // ' atom/b-cm'
write(unit_,*) ' Atom Density = ' // trim(to_str(mat % density)) &
// ' atom/b-cm'
! Write atom density for each nuclide in material
write(ou,*) ' Nuclides:'
write(unit_,*) ' Nuclides:'
do i = 1, mat % n_nuclides
nuc => nuclides(mat % nuclide(i))
density = mat % atom_density(i)
string = ' ' // trim(nuc % name) // ' = ' // &
& trim(real_to_str(density)) // ' atom/b-cm'
write(ou,*) trim(string)
trim(to_str(density)) // ' atom/b-cm'
write(unit_,*) trim(string)
end do
! Write information on S(a,b) table
if (mat % has_sab_table) then
write(ou,*) ' S(a,b) table = ' // trim(mat % sab_name)
write(unit_,*) ' S(a,b) table = ' // trim(mat % sab_name)
end if
write(ou,*)
write(unit_,*)
end subroutine print_material
@ -501,29 +548,37 @@ contains
! PRINT_TALLY displays the attributes of a tally
!===============================================================================
subroutine print_tally(t)
subroutine print_tally(t, unit)
type(TallyObject), pointer :: t
integer, optional :: unit
integer :: i
integer :: id
character(MAX_LINE_LEN) :: string
integer :: i
integer :: id
integer :: unit_
character(MAX_LINE_LEN) :: string
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Universe), pointer :: u => null()
type(Material), pointer :: m => null()
type(StructuredMesh), pointer :: sm => null()
write(ou,*) 'Tally ' // int_to_str(t % id)
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
write(unit_,*) 'Tally ' // to_str(t % id)
if (t % n_bins(T_CELL) > 0) then
string = ""
do i = 1, t % n_bins(T_CELL)
id = t % cell_bins(i) % scalar
c => cells(id)
string = trim(string) // ' ' // trim(int_to_str(c % id))
string = trim(string) // ' ' // trim(to_str(c % id))
end do
write(ou, *) ' Cell Bins:' // trim(string)
write(unit_, *) ' Cell Bins:' // trim(string)
end if
if (t % n_bins(T_SURFACE) > 0) then
@ -531,9 +586,9 @@ contains
do i = 1, t % n_bins(T_SURFACE)
id = t % surface_bins(i) % scalar
s => surfaces(id)
string = trim(string) // ' ' // trim(int_to_str(s % id))
string = trim(string) // ' ' // trim(to_str(s % id))
end do
write(ou, *) ' Surface Bins:' // trim(string)
write(unit_, *) ' Surface Bins:' // trim(string)
end if
if (t % n_bins(T_UNIVERSE) > 0) then
@ -541,9 +596,9 @@ contains
do i = 1, t % n_bins(T_UNIVERSE)
id = t % universe_bins(i) % scalar
u => universes(id)
string = trim(string) // ' ' // trim(int_to_str(u % id))
string = trim(string) // ' ' // trim(to_str(u % id))
end do
write(ou, *) ' Material Bins:' // trim(string)
write(unit_, *) ' Material Bins:' // trim(string)
end if
if (t % n_bins(T_MATERIAL) > 0) then
@ -551,20 +606,20 @@ contains
do i = 1, t % n_bins(T_MATERIAL)
id = t % material_bins(i) % scalar
m => materials(id)
string = trim(string) // ' ' // trim(int_to_str(m % id))
string = trim(string) // ' ' // trim(to_str(m % id))
end do
write(ou, *) ' Material Bins:' // trim(string)
write(unit_, *) ' Material Bins:' // trim(string)
end if
if (t % n_bins(T_MESH) > 0) then
string = ""
id = t % mesh
sm => meshes(id)
string = trim(string) // ' ' // trim(int_to_str(sm % dimension(1)))
string = trim(string) // ' ' // trim(to_str(sm % dimension(1)))
do i = 2, sm % n_dimension
string = trim(string) // ' x ' // trim(int_to_str(sm % dimension(i)))
string = trim(string) // ' x ' // trim(to_str(sm % dimension(i)))
end do
write(ou, *) ' Mesh Bins:' // trim(string)
write(unit_, *) ' Mesh Bins:' // trim(string)
end if
if (t % n_bins(T_CELLBORN) > 0) then
@ -572,27 +627,27 @@ contains
do i = 1, t % n_bins(T_CELLBORN)
id = t % cellborn_bins(i) % scalar
c => cells(id)
string = trim(string) // ' ' // trim(int_to_str(c % id))
string = trim(string) // ' ' // trim(to_str(c % id))
end do
write(ou, *) ' Birth Region Bins:' // trim(string)
write(unit_, *) ' Birth Region Bins:' // trim(string)
end if
if (t % n_bins(T_ENERGYIN) > 0) then
string = ""
do i = 1, t % n_bins(T_ENERGYIN) + 1
string = trim(string) // ' ' // trim(real_to_str(&
string = trim(string) // ' ' // trim(to_str(&
t % energy_in(i)))
end do
write(ou,*) ' Incoming Energy Bins:' // trim(string)
write(unit_,*) ' Incoming Energy Bins:' // trim(string)
end if
if (t % n_bins(T_ENERGYOUT) > 0) then
string = ""
do i = 1, t % n_bins(T_ENERGYOUT) + 1
string = trim(string) // ' ' // trim(real_to_str(&
string = trim(string) // ' ' // trim(to_str(&
t % energy_out(i)))
end do
write(ou,*) ' Outgoing Energy Bins:' // trim(string)
write(unit_,*) ' Outgoing Energy Bins:' // trim(string)
end if
if (t % n_macro_bins > 0) then
@ -613,9 +668,9 @@ contains
string = trim(string) // ' nu-fission'
end select
end do
write(ou,*) ' Macro Reactions:' // trim(string)
write(unit_,*) ' Macro Reactions:' // trim(string)
end if
write(ou,*)
write(unit_,*)
end subroutine print_tally
@ -632,37 +687,129 @@ contains
type(Universe), pointer :: u => null()
type(Lattice), pointer :: l => null()
! print summary of surfaces
call header("SURFACE SUMMARY", unit=UNIT_SUMMARY)
do i = 1, n_surfaces
s => surfaces(i)
call print_surface(s, unit=UNIT_SUMMARY)
end do
! print summary of cells
call header("CELL SUMMARY")
call header("CELL SUMMARY", unit=UNIT_SUMMARY)
do i = 1, n_cells
c => cells(i)
call print_cell(c)
call print_cell(c, unit=UNIT_SUMMARY)
end do
! print summary of universes
call header("UNIVERSE SUMMARY")
call header("UNIVERSE SUMMARY", unit=UNIT_SUMMARY)
do i = 1, n_universes
u => universes(i)
call print_universe(u)
call print_universe(u, unit=UNIT_SUMMARY)
end do
! print summary of lattices
if (n_lattices > 0) then
call header("LATTICE SUMMARY")
call header("LATTICE SUMMARY", unit=UNIT_SUMMARY)
do i = 1, n_lattices
l => lattices(i)
call print_lattice(l)
call print_lattice(l, unit=UNIT_SUMMARY)
end do
end if
! print summary of surfaces
call header("SURFACE SUMMARY")
do i = 1, n_surfaces
s => surfaces(i)
call print_surface(s)
end subroutine print_geometry
!===============================================================================
! PRINT_NUCLIDE displays information about a continuous-energy neutron
! cross_section table and its reactions and secondary angle/energy distributions
!===============================================================================
subroutine print_nuclide(nuc, unit)
type(Nuclide), pointer :: nuc
integer, optional :: unit
integer :: i
integer :: unit_
integer :: size_total
integer :: size_xs
integer :: size_angle
integer :: size_energy
type(Reaction), pointer :: rxn => null()
type(UrrData), pointer :: urr => null()
! set default unit for writing information
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Determine size of cross-sections
size_xs = (5 + nuc % n_reaction) * nuc % n_grid * 8
size_total = size_xs
! Basic nuclide information
write(unit_,*) 'Nuclide ' // trim(nuc % name)
write(unit_,*) ' zaid = ' // trim(to_str(nuc % zaid))
write(unit_,*) ' awr = ' // trim(to_str(nuc % awr))
write(unit_,*) ' kT = ' // trim(to_str(nuc % kT))
write(unit_,*) ' # of grid points = ' // trim(to_str(nuc % n_grid))
write(unit_,*) ' Fissionable = ', nuc % fissionable
write(unit_,*) ' # of fission reactions = ' // trim(to_str(nuc % n_fission))
write(unit_,*) ' # of reactions = ' // trim(to_str(nuc % n_reaction))
write(unit_,*) ' Size of cross sections = ' // trim(to_str(&
size_xs)) // ' bytes'
write(unit_,*) ' Reaction Q-value Mult IE size(angle) size(energy)'
do i = 1, nuc % n_reaction
! Information on each reaction
rxn => nuc % reactions(i)
! Determine size of angle distribution
if (rxn % has_angle_dist) then
size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8
else
size_angle = 0
end if
! Determine size of energy distribution
if (rxn % has_energy_dist) then
size_energy = size(rxn % edist % data) * 8
else
size_energy = 0
end if
write(unit_,'(3X,A11,1X,F8.3,2X,I4,2X,I6,1X,I11,1X,I11)') &
reaction_name(rxn % MT), rxn % Q_value, rxn % TY, rxn % IE, &
size_angle, size_energy
! Accumulate data size
size_total = size_total + size_angle + size_energy
end do
end subroutine print_geometry
! Write information about URR probability tables
if (nuc % urr_present) then
urr => nuc % urr_data
write(unit_,*) ' Unresolved resonance probability table:'
write(unit_,*) ' # of energies = ' // trim(to_str(urr % n_energy))
write(unit_,*) ' # of probabilities = ' // trim(to_str(urr % n_prob))
write(unit_,*) ' Interpolation = ' // trim(to_str(urr % interp))
write(unit_,*) ' Inelastic flag = ' // trim(to_str(urr % inelastic_flag))
write(unit_,*) ' Absorption flag = ' // trim(to_str(urr % absorption_flag))
write(unit_,*) ' Multiply by smooth? ', urr % multiply_smooth
write(unit_,*) ' Min energy = ', trim(to_str(urr % energy(1)))
write(unit_,*) ' Max energy = ', trim(to_str(urr % energy(urr % n_energy)))
end if
! Write total memory used
write(unit_,*) ' Total memory used = ' // trim(to_str(size_total)) &
// ' bytes'
! Blank line at end of nuclide
write(unit_,*)
end subroutine print_nuclide
!===============================================================================
! PRINT_SUMMARY displays summary information about the problem about to be run
@ -677,65 +824,63 @@ contains
type(TallyObject), pointer :: t => null()
! Display problem summary
call header("PROBLEM SUMMARY")
call header("PROBLEM SUMMARY", unit=UNIT_SUMMARY)
if (problem_type == PROB_CRITICALITY) then
write(ou,100) 'Problem type:', 'Criticality'
write(ou,101) 'Number of Cycles:', n_cycles
write(ou,101) 'Number of Inactive Cycles:', n_inactive
write(UNIT_SUMMARY,100) 'Problem type:', 'Criticality'
write(UNIT_SUMMARY,101) 'Number of Cycles:', n_cycles
write(UNIT_SUMMARY,101) 'Number of Inactive Cycles:', n_inactive
elseif (problem_type == PROB_SOURCE) then
write(ou,100) 'Problem type:', 'External Source'
write(UNIT_SUMMARY,100) 'Problem type:', 'External Source'
end if
write(ou,101) 'Number of Particles:', n_particles
write(UNIT_SUMMARY,101) 'Number of Particles:', n_particles
! Display geometry summary
call header("GEOMETRY SUMMARY")
write(ou,101) 'Number of Cells:', n_cells
write(ou,101) 'Number of Surfaces:', n_surfaces
write(ou,101) 'Number of Materials:', n_materials
call header("GEOMETRY SUMMARY", unit=UNIT_SUMMARY)
write(UNIT_SUMMARY,101) 'Number of Cells:', n_cells
write(UNIT_SUMMARY,101) 'Number of Surfaces:', n_surfaces
write(UNIT_SUMMARY,101) 'Number of Materials:', n_materials
! print summary of all geometry
call print_geometry()
! print summary of materials
call header("MATERIAL SUMMARY")
call header("MATERIAL SUMMARY", unit=UNIT_SUMMARY)
do i = 1, n_materials
m => materials(i)
call print_material(m)
call print_material(m, unit=UNIT_SUMMARY)
end do
! print summary of tallies
if (n_tallies > 0) then
call header("TALLY SUMMARY")
call header("TALLY SUMMARY", unit=UNIT_SUMMARY)
do i = 1, n_tallies
t=> tallies(i)
call print_tally(t)
call print_tally(t, unit=UNIT_SUMMARY)
end do
end if
! print summary of unionized energy grid
call header("UNIONIZED ENERGY GRID")
write(ou,*) "Points on energy grid: " // trim(int_to_str(n_grid))
write(ou,*) "Extra storage required: " // trim(int_to_str(&
call header("UNIONIZED ENERGY GRID", unit=UNIT_SUMMARY)
write(UNIT_SUMMARY,*) "Points on energy grid: " // trim(to_str(n_grid))
write(UNIT_SUMMARY,*) "Extra storage required: " // trim(to_str(&
n_grid*n_nuclides_total*4)) // " bytes"
! print summary of variance reduction
call header("VARIANCE REDUCTION")
call header("VARIANCE REDUCTION", unit=UNIT_SUMMARY)
if (survival_biasing) then
write(ou,100) "Survival Biasing:", "on"
write(UNIT_SUMMARY,100) "Survival Biasing:", "on"
else
write(ou,100) "Survival Biasing:", "off"
write(UNIT_SUMMARY,100) "Survival Biasing:", "off"
end if
string = real_to_str(weight_cutoff)
write(ou,100) "Weight Cutoff:", trim(string)
string = real_to_str(weight_survive)
write(ou,100) "Survival weight:", trim(string)
write(ou,*)
string = to_str(weight_cutoff)
write(UNIT_SUMMARY,100) "Weight Cutoff:", trim(string)
string = to_str(weight_survive)
write(UNIT_SUMMARY,100) "Survival weight:", trim(string)
! Format descriptor for columns
100 format (1X,A,T35,A)
101 format (1X,A,T35,I11)
end subroutine print_summary
!===============================================================================
@ -748,16 +893,16 @@ contains
call header("PLOTTING SUMMARY")
! Print plotting origin
write(ou,100) "Plotting Origin:", trim(real_to_str(plot_origin(1))) // &
" " // trim(real_to_str(plot_origin(2))) // " " // &
trim(real_to_str(plot_origin(3)))
write(ou,100) "Plotting Origin:", trim(to_str(plot_origin(1))) // &
" " // trim(to_str(plot_origin(2))) // " " // &
trim(to_str(plot_origin(3)))
! Print plotting width
write(ou,100) "Plotting Width:", trim(real_to_str(plot_width(1))) // &
" " // trim(real_to_str(plot_width(2)))
write(ou,100) "Plotting Width:", trim(to_str(plot_width(1))) // &
" " // trim(to_str(plot_width(2)))
! Print pixel width
write(ou,100) "Pixel Width:", trim(real_to_str(pixel))
write(ou,100) "Pixel Width:", trim(to_str(pixel))
write(ou,*)
! Format descriptor for columns
@ -799,7 +944,7 @@ contains
! display calculate rate and final keff
total_particles = n_particles * n_cycles
speed = real(total_particles) / time_compute % elapsed
string = real_to_str(speed)
string = to_str(speed)
write(ou,101) "Calculation Rate", trim(string)
write(ou,102) "Final Keff", keff, keff_std
write(ou,*)
@ -811,4 +956,30 @@ contains
end subroutine print_runtime
!===============================================================================
! CREATE_SUMMARY_FILE opens the summary.out file for logging information about
! the simulation
!===============================================================================
subroutine create_summary_file()
integer :: io_error
logical :: file_exists ! does log file already exist?
character(MAX_FILE_LEN) :: path ! path of summary file
! Create filename for log file
path = "summary.out"
! Check if log file already exists
inquire(FILE=path, EXIST=file_exists)
if (file_exists) then
! Possibly copy old log file
end if
! Open log file for writing
open(UNIT=UNIT_SUMMARY, FILE=path, STATUS='replace', &
ACTION='write', IOSTAT=io_error)
end subroutine create_summary_file
end module output

361
src/physics.F90
View file

@ -1,21 +1,22 @@
module physics
use ace_header, only: Nuclide, Reaction, DistEnergy
use constants
use cross_section_header, only: Nuclide, Reaction, DistEnergy
use endf, only: reaction_name, is_fission, is_scatter
use error, only: fatal_error, warning
use fission, only: nu_total, nu_prompt, nu_delayed
use geometry, only: find_cell, distance_to_boundary, &
cross_surface, cross_lattice
use geometry_header, only: Universe, BASE_UNIVERSE
use cross_section, only: calculate_xs
use endf, only: reaction_name, is_fission, is_scatter
use error, only: fatal_error, warning
use fission, only: nu_total, nu_prompt, nu_delayed
use geometry, only: find_cell, distance_to_boundary, cross_surface, &
cross_lattice
use geometry_header, only: Universe, BASE_UNIVERSE
use global
use interpolation, only: interpolate_tab1
use output, only: write_message
use particle_header, only: Particle, LocalCoord
use random_lcg, only: prn
use search, only: binary_search
use string, only: int_to_str, real_to_str
use tally, only: score_tally, score_surface_current
use interpolation, only: interpolate_tab1
use output, only: write_message
use particle_header, only: Particle, LocalCoord
use random_lcg, only: prn
use search, only: binary_search
use string, only: to_str
use tally, only: score_tally, score_surface_current
implicit none
@ -29,22 +30,21 @@ contains
type(Particle), pointer :: p
integer :: surface_crossed ! surface which particle is on
integer :: last_cell ! most recent cell particle was in
integer :: n_event ! number of collisions/crossings
real(8) :: d_boundary ! distance to nearest boundary
real(8) :: d_collision ! sampled distance to collision
real(8) :: distance ! distance particle travels
logical :: found_cell ! found cell which particle is in?
logical :: lattice_crossed ! is surface crossing in lattice?
integer :: surface_crossed ! surface which particle is on
integer :: lattice_crossed ! lattice boundary which particle crossed
integer :: last_cell ! most recent cell particle was in
integer :: n_event ! number of collisions/crossings
real(8) :: d_boundary ! distance to nearest boundary
real(8) :: d_collision ! sampled distance to collision
real(8) :: distance ! distance particle travels
logical :: found_cell ! found cell which particle is in?
type(LocalCoord), pointer :: coord => null()
if (p % coord % cell == NONE) then
call find_cell(p, found_cell)
! Particle couldn't be located
if (.not. found_cell) then
write(message, '(A,3ES11.3)') &
"Could not locate cell for particle at: ", p % coord0 % xyz
message = "Could not locate particle " // trim(to_str(p % id))
call fatal_error()
end if
@ -53,12 +53,12 @@ contains
end if
if (verbosity >= 9 .or. trace) then
message = "Simulating Particle " // trim(int_to_str(p % id))
message = "Simulating Particle " // trim(to_str(p % id))
call write_message()
end if
if (verbosity >= 10 .or. trace) then
message = " Born in cell " // trim(int_to_str(&
message = " Born in cell " // trim(to_str(&
cells(p % coord % cell) % id))
call write_message()
end if
@ -94,9 +94,9 @@ contains
if (d_collision > d_boundary) then
last_cell = p % coord % cell
p % coord % cell = NONE
if (lattice_crossed) then
if (lattice_crossed /= NONE) then
p % surface = NONE
call cross_lattice(p)
call cross_lattice(p, lattice_crossed)
else
p % surface = surface_crossed
call cross_surface(p, last_cell)
@ -125,8 +125,8 @@ contains
! If particle has too many events, display warning and kill it
n_event = n_event + 1
if (n_event == MAX_EVENTS) then
message = "Particle " // trim(int_to_str(p%id)) // " underwent " &
// "maximum number of events."
message = "Particle " // trim(to_str(p%id)) // " underwent maximum &
&number of events."
call warning()
p % alive = .false.
end if
@ -135,259 +135,6 @@ contains
end subroutine transport
!===============================================================================
! CALCULATE_XS determines the macroscopic cross sections for the material the
! particle is currently traveling through.
!===============================================================================
subroutine calculate_xs(p)
type(Particle), pointer :: p
integer :: i ! loop index over nuclides
integer :: index_nuclide ! index into nuclides array
integer :: index_sab ! index into sab_tables array
real(8) :: atom_density ! atom density of a nuclide
real(8) :: sab_threshold ! threshold for S(a,b) table
type(Material), pointer :: mat => null() ! current material
! Set all material macroscopic cross sections to zero
material_xs % total = ZERO
material_xs % elastic = ZERO
material_xs % absorption = ZERO
material_xs % fission = ZERO
material_xs % nu_fission = ZERO
mat => materials(p % material)
! Find energy index on unionized grid
call find_energy_index(p)
! Check if there's an S(a,b) table for this material
if (mat % has_sab_table) then
sab_threshold = sab_tables(mat % sab_table) % threshold_inelastic
else
sab_threshold = ZERO
end if
! Add contribution from each nuclide in material
do i = 1, mat % n_nuclides
! Determine microscopic cross sections for this nuclide
index_nuclide = mat % nuclide(i)
! Determine whether to use S(a,b) based on energy of particle and whether
! the nuclide matches the S(a,b) table
if (p % E < sab_threshold .and. i == mat % sab_nuclide) then
index_sab = mat % sab_table
else
index_sab = 0
end if
! Calculate microscopic cross section for this nuclide
if (p % E /= micro_xs(index_nuclide) % last_E) then
call calculate_nuclide_xs(p, index_nuclide, index_sab)
end if
! Copy atom density of nuclide in material
atom_density = mat % atom_density(i)
! Add contributions to material macroscopic total cross section
material_xs % total = material_xs % total + &
atom_density * micro_xs(index_nuclide) % total
! Add contributions to material macroscopic scattering cross section
material_xs % elastic = material_xs % elastic + &
atom_density * micro_xs(index_nuclide) % elastic
! Add contributions to material macroscopic absorption cross section
material_xs % absorption = material_xs % absorption + &
atom_density * micro_xs(index_nuclide) % absorption
! Add contributions to material macroscopic fission cross section
material_xs % fission = material_xs % fission + &
atom_density * micro_xs(index_nuclide) % fission
! Add contributions to material macroscopic nu-fission cross section
material_xs % nu_fission = material_xs % nu_fission + &
atom_density * micro_xs(index_nuclide) % nu_fission
end do
end subroutine calculate_xs
!===============================================================================
! CALCULATE_NUCLIDE_XS determines microscopic cross sections for a nuclide of a
! given index in the nuclides array at the energy of the given particle
!===============================================================================
subroutine calculate_nuclide_xs(p, index_nuclide, index_sab)
type(Particle), pointer :: p
integer, intent(in) :: index_nuclide ! index into nuclides array
integer, intent(in) :: index_sab ! index into sab_tables array
integer :: IE ! index on nuclide energy grid
integer :: IE_sab ! index on S(a,b) energy grid
real(8) :: f ! interp factor on nuclide energy grid
real(8) :: f_sab ! interp factor on S(a,b) energy grid
real(8) :: inelastic ! S(a,b) inelastic cross section
real(8) :: elastic ! S(a,b) elastic cross section
type(Nuclide), pointer :: nuc => null()
type(SAB_Table), pointer :: sab => null()
! Set pointer to nuclide
nuc => nuclides(index_nuclide)
! TODO: If not using unionized energy grid, we need to find the index on the
! nuclide energy grid using lethargy mapping or whatever other technique
! search nuclide energy grid
IE = nuc % grid_index(p % IE)
f = (p%E - nuc%energy(IE))/(nuc%energy(IE+1) - nuc%energy(IE))
micro_xs(index_nuclide) % index_grid = IE
micro_xs(index_nuclide) % interp_factor = f
! Initialize sab treatment to false
micro_xs(index_nuclide) % use_sab = .false.
micro_xs(index_nuclide) % elastic_sab = ZERO
! Initialize nuclide cross-sections to zero
micro_xs(index_nuclide) % fission = ZERO
micro_xs(index_nuclide) % nu_fission = ZERO
! Calculate microscopic nuclide total cross section
micro_xs(index_nuclide) % total = &
(ONE-f) * nuc % total(IE) + f * nuc % total(IE+1)
! Calculate microscopic nuclide total cross section
micro_xs(index_nuclide) % elastic = &
(ONE-f) * nuc % elastic(IE) + f * nuc % elastic(IE+1)
! Calculate microscopic nuclide absorption cross section
micro_xs(index_nuclide) % absorption = &
(ONE-f) * nuc % absorption(IE) + f * nuc % absorption(IE+1)
if (nuc % fissionable) then
! Calculate microscopic nuclide total cross section
micro_xs(index_nuclide) % fission = &
(ONE-f) * nuc % fission(IE) + f * nuc % fission(IE+1)
! Calculate microscopic nuclide nu-fission cross section
micro_xs(index_nuclide) % nu_fission = &
(ONE-f) * nuc % nu_fission(IE) + f * nuc % nu_fission(IE+1)
end if
! If there is S(a,b) data for this nuclide, we need to do a few
! things. Since the total cross section was based on non-S(a,b) data, we
! need to correct it by subtracting the non-S(a,b) elastic cross section and
! then add back in the calculated S(a,b) elastic+inelastic cross section.
if (index_sab > 0) then
micro_xs(index_nuclide) % use_sab = .true.
! Get pointer to S(a,b) table
sab => sab_tables(index_sab)
! Get index and interpolation factor for inelastic grid
if (p%E < sab % inelastic_e_in(1)) then
IE_sab = 1
f_sab = ZERO
else
IE_sab = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, p%E)
f_sab = (p%E - sab%inelastic_e_in(IE_sab)) / &
(sab%inelastic_e_in(IE_sab+1) - sab%inelastic_e_in(IE_sab))
end if
! Calculate S(a,b) inelastic scattering cross section
inelastic = (ONE-f_sab) * sab % inelastic_sigma(IE_sab) + f_sab * &
sab % inelastic_sigma(IE_sab + 1)
! Check for elastic data
if (p % E < sab % threshold_elastic) then
! Determine whether elastic scattering is given in the coherent or
! incoherent approximation. For coherent, the cross section is
! represented as P/E whereas for incoherent, it is simply P
if (sab % elastic_mode == SAB_ELASTIC_EXACT) then
if (p % E < sab % elastic_e_in(1)) then
! If energy is below that of the lowest Bragg peak, the elastic
! cross section will be zero
elastic = ZERO
else
IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E)
elastic = sab % elastic_P(IE_sab) / p % E
end if
else
! Determine index on elastic energy grid
if (p % E < sab % elastic_e_in(1)) then
IE_sab = 1
else
IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E)
end if
! Get interpolation factor for elastic grid
f_sab = (p%E - sab%elastic_e_in(IE_sab))/(sab%elastic_e_in(IE_sab+1) - &
sab%elastic_e_in(IE_sab))
! Calculate S(a,b) elastic scattering cross section
elastic = (ONE-f_sab) * sab % elastic_P(IE_sab) + f_sab * &
sab % elastic_P(IE_sab + 1)
end if
else
! No elastic data
elastic = ZERO
end if
! Correct total and elastic cross sections
micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % total - &
micro_xs(index_nuclide) % elastic + inelastic + elastic
micro_xs(index_nuclide) % elastic = inelastic + elastic
! Store S(a,b) elastic cross section for sampling later
micro_xs(index_nuclide) % elastic_sab = elastic
end if
! Set last evaluated energy
micro_xs(index_nuclide) % last_E = p % E
end subroutine calculate_nuclide_xs
!===============================================================================
! FIND_ENERGY_INDEX determines the index on the union energy grid and the
! interpolation factor for a particle at a certain energy
!===============================================================================
subroutine find_energy_index(p)
type(Particle), pointer :: p
integer :: IE ! index on union energy grid
real(8) :: E ! energy of particle
real(8) :: interp ! interpolation factor
! copy particle's energy
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 union energy grid.
if (E < e_grid(1)) then
IE = 1
elseif (E > e_grid(n_grid)) then
IE = n_grid - 1
else
IE = binary_search(e_grid, n_grid, E)
end if
! calculate the interpolation factor -- note this will be outside of [0,1)
! for a particle outside the energy range of the union grid
interp = (E - e_grid(IE))/(e_grid(IE+1) - e_grid(IE))
! set particle attributes
p % IE = IE
p % interp = interp
end subroutine find_energy_index
!===============================================================================
! COLLISION samples a nuclide and reaction and then calls the appropriate
! routine for that reaction
@ -420,7 +167,7 @@ contains
! Display information about collision
if (verbosity >= 10 .or. trace) then
message = " " // trim(reaction_name(MT)) // ". Energy = " // &
trim(real_to_str(p % E * 1e6_8)) // " eV."
trim(to_str(p % E * 1e6_8)) // " eV."
call write_message()
end if
@ -450,9 +197,6 @@ contains
! Reset number of particles banked during collision
p % n_bank = 0
! find energy index, interpolation factor
call find_energy_index(p)
end subroutine collision
!===============================================================================
@ -571,7 +315,7 @@ contains
! add to cumulative probability
prob = prob + ((ONE-f)*rxn%sigma(IE-rxn%IE+1) &
& + f*(rxn%sigma(IE-rxn%IE+2)))
+ f*(rxn%sigma(IE-rxn%IE+2)))
end do
else
! For nuclides with only total fission reaction, get a pointer to
@ -596,7 +340,7 @@ contains
! add to cumulative probability
if (nuc % has_partial_fission) then
prob = prob + ((ONE-f)*rxn%sigma(IE-rxn%IE+1) &
& + f*(rxn%sigma(IE-rxn%IE+2)))
+ f*(rxn%sigma(IE-rxn%IE+2)))
else
prob = prob + micro_xs(index_nuclide) % fission
end if
@ -672,7 +416,7 @@ contains
if (i > nuc % n_reaction) then
message = "Did not sample any reaction for nuclide " // &
trim(nuc % name) // " on material " // &
trim(int_to_str(mat % id))
trim(to_str(mat % id))
call fatal_error()
end if
@ -692,7 +436,7 @@ contains
! add to cumulative probability
prob = prob + ((ONE-f)*rxn%sigma(IE-rxn%IE+1) &
& + f*(rxn%sigma(IE-rxn%IE+2)))
+ f*(rxn%sigma(IE-rxn%IE+2)))
end do
! Perform collision physics for inelastics scattering
@ -1023,7 +767,7 @@ contains
! determine velocity vector of target nucleus
v_target(1) = u*vt
v_target(2) = v*vt
v_target(3) = W*vt
v_target(3) = w*vt
end subroutine sample_target_velocity
@ -1262,7 +1006,7 @@ contains
! determine outgoing energy in lab
E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) &
& / ((A+ONE)*(A+ONE))
/ ((A+ONE)*(A+ONE))
! determine outgoing angle in lab
mu = mu * sqrt(E_cm/E) + ONE/(A+ONE) * sqrt(E_in/E)
@ -1338,7 +1082,7 @@ contains
else
i = binary_search(rxn % adist % energy, n, E)
r = (E - rxn % adist % energy(i)) / &
& (rxn % adist % energy(i+1) - rxn % adist % energy(i))
(rxn % adist % energy(i+1) - rxn % adist % energy(i))
end if
! Sample between the ith and (i+1)th bin
@ -1380,8 +1124,7 @@ contains
mu = mu0 + (xi - c_k)/p0
elseif (interp == LINEAR_LINEAR) then
! Linear-linear interpolation -- not sure how you come about the
! formula given in the MCNP manual
! Linear-linear interpolation
p1 = rxn % adist % data(lc + NP + k+1)
mu1 = rxn % adist % data(lc + k+1)
@ -1392,7 +1135,7 @@ contains
mu = mu0 + (sqrt(p0*p0 + 2*frac*(xi - c_k))-p0)/frac
end if
else
message = "Unknown interpolation type: " // trim(int_to_str(interp))
message = "Unknown interpolation type: " // trim(to_str(interp))
call fatal_error()
end if
@ -1404,7 +1147,7 @@ contains
end if
else
message = "Unknown angular distribution type: " // trim(int_to_str(type))
message = "Unknown angular distribution type: " // trim(to_str(type))
call fatal_error()
end if
@ -1559,7 +1302,7 @@ contains
lc = 2 + 2*NR
i = binary_search(edist % data(lc+1), NE, E_in)
r = (E_in - edist%data(lc+i)) / &
& (edist%data(lc+i+1) - edist%data(lc+i))
(edist%data(lc+i+1) - edist%data(lc+i))
! Sample outgoing energy bin
r1 = prn()
@ -1635,7 +1378,7 @@ contains
else
i = binary_search(edist % data(lc+1), NE, E_in)
r = (E_in - edist%data(lc+i)) / &
& (edist%data(lc+i+1) - edist%data(lc+i))
(edist%data(lc+i+1) - edist%data(lc+i))
end if
! Sample between the ith and (i+1)th bin
@ -1708,10 +1451,10 @@ contains
E_out = E_l_k + (r1 - c_k)/p_l_k
else
E_out = E_l_k + (sqrt(p_l_k*p_l_k + 2*frac*(r1 - c_k)) - &
& p_l_k)/frac
p_l_k)/frac
end if
else
message = "Unknown interpolation type: " // trim(int_to_str(INTT))
message = "Unknown interpolation type: " // trim(to_str(INTT))
call fatal_error()
end if
@ -1861,7 +1604,7 @@ contains
else
i = binary_search(edist % data(lc+1), NE, E_in)
r = (E_in - edist%data(lc+i)) / &
& (edist%data(lc+i+1) - edist%data(lc+i))
(edist%data(lc+i+1) - edist%data(lc+i))
end if
! Sample between the ith and (i+1)th bin
@ -1940,7 +1683,7 @@ contains
E_out = E_l_k + (r1 - c_k)/p_l_k
else
E_out = E_l_k + (sqrt(p_l_k*p_l_k + 2*frac*(r1 - c_k)) - &
& p_l_k)/frac
p_l_k)/frac
end if
! Determine Kalbach-Mann parameters
@ -1952,7 +1695,7 @@ contains
KM_R = R_k + (R_k1 - R_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k)
KM_A = A_k + (A_k1 - A_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k)
else
message = "Unknown interpolation type: " // trim(int_to_str(INTT))
message = "Unknown interpolation type: " // trim(to_str(INTT))
call fatal_error()
end if
@ -2004,7 +1747,7 @@ contains
else
i = binary_search(edist % data(lc+1), NE, E_in)
r = (E_in - edist%data(lc+i)) / &
& (edist%data(lc+i+1) - edist%data(lc+i))
(edist%data(lc+i+1) - edist%data(lc+i))
end if
! Sample between the ith and (i+1)th bin
@ -2079,10 +1822,10 @@ contains
E_out = E_l_k + (r1 - c_k)/p_l_k
else
E_out = E_l_k + (sqrt(p_l_k*p_l_k + 2*frac*(r1 - c_k)) - &
& p_l_k)/frac
p_l_k)/frac
end if
else
message = "Unknown interpolation type: " // trim(int_to_str(INTT))
message = "Unknown interpolation type: " // trim(to_str(INTT))
call fatal_error()
end if
@ -2135,7 +1878,7 @@ contains
mu_out = mu_k + (sqrt(p_k*p_k + 2*frac*(r3 - c_k))-p_k)/frac
end if
else
message = "Unknown interpolation type: " // trim(int_to_str(JJ))
message = "Unknown interpolation type: " // trim(to_str(JJ))
call fatal_error()
end if

6
src/plot.F90
View file

@ -23,6 +23,7 @@ contains
integer :: i ! loop index
integer :: surface_crossed ! surface which particle is on
integer :: lattice_crossed ! is surface crossing in lattice?
integer :: last_cell ! most recent cell particle was in
integer :: enter_surface ! entrance surface
real(8) :: xyz(3) ! starting coordinates
@ -31,7 +32,6 @@ contains
real(8) :: d ! distance to boundary
real(8) :: distance ! distance particle travels
logical :: found_cell ! found cell which particle is in?
logical :: lattice_crossed ! is surface crossing in lattice?
character(MAX_LINE_LEN) :: path_plot ! unit for binary plot file
type(Cell), pointer :: c => null()
type(Universe), pointer :: univ => null()
@ -161,9 +161,9 @@ contains
write(UNIT=UNIT_PLOT) p % coord0 % xyz, last_cell
p % coord % cell = 0
if (lattice_crossed) then
if (lattice_crossed /= NONE) then
p % surface = NONE
call cross_lattice(p)
call cross_lattice(p, lattice_crossed)
else
p % surface = surface_crossed
call cross_surface(p, last_cell)

2
src/random_lcg.F90
View file

@ -23,7 +23,7 @@ module random_lcg
contains
!===============================================================================
! PRN generates a pseudo-random number using a linean congruential generator
! PRN generates a pseudo-random number using a linear congruential generator
!===============================================================================
function prn() result(pseudo_rn)

26
src/source.F90
View file

@ -1,15 +1,15 @@
module source
use bank_header, only: Bank
use constants, only: ONE, MAX_LINE_LEN
use cross_section_header, only: Nuclide
use error, only: fatal_error
use ace_header, only: Nuclide
use bank_header, only: Bank
use constants, only: ONE, MAX_LINE_LEN
use error, only: fatal_error
use global
use output, only: write_message
use particle_header, only: Particle, initialize_particle
use physics, only: watt_spectrum
use random_lcg, only: prn, set_particle_seed
use string, only: int_to_str
use output, only: write_message
use particle_header, only: Particle, initialize_particle
use physics, only: watt_spectrum
use random_lcg, only: prn, set_particle_seed
use string, only: to_str
implicit none
@ -51,7 +51,7 @@ contains
bytes = maxwork * 64 / 8
!#endif
message = "Could not allocate source bank. Attempted to allocate " &
// trim(int_to_str(bytes)) // " bytes."
// trim(to_str(bytes)) // " bytes."
call fatal_error()
end if
@ -64,7 +64,7 @@ contains
bytes = 3 * maxwork * 64 / 8
!#endif
message = "Could not allocate fission bank. Attempted to allocate " &
// trim(int_to_str(bytes)) // " bytes."
// trim(to_str(bytes)) // " bytes."
call fatal_error()
end if
@ -72,6 +72,10 @@ contains
if (external_source%type == SRC_BOX) then
p_min = external_source%values(1:3)
p_max = external_source%values(4:6)
else
message = "Unsupported external source type: " // &
to_str(external_source%type)
call fatal_error()
end if
! Initialize first cycle source bank

6
src/string.F90
View file

@ -6,8 +6,8 @@ module string
implicit none
interface int_to_str
module procedure int4_to_str, int8_to_str
interface to_str
module procedure int4_to_str, int8_to_str, real_to_str
end interface
contains
@ -348,7 +348,7 @@ contains
! Determine number of digits to right of decimal point
index_decimal = index(string, '.')
index_exponent = max(index(string, 'd'), index(string, 'D'), &
& index(string, 'e'), index(string, 'E'))
index(string, 'e'), index(string, 'E'))
if (index_decimal > 0) then
if (index_exponent > 0) then
d = index_exponent - index_decimal - 1

165
src/tally.F90
View file

@ -1,20 +1,18 @@
module tally
use ISO_FORTRAN_ENV
use constants
use error, only: fatal_error
use global
use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices
use mesh_header, only: StructuredMesh
use output, only: write_message, header
use output, only: header
use search, only: binary_search
use string, only: int_to_str, real_to_str
use string, only: to_str
use tally_header, only: TallyScore, TallyMapItem, TallyMapElement
#ifdef MPI
use intercycle, only: reduce_tallies
use mpi
use mpi_routines, only: reduce_tallies
#endif
implicit none
@ -23,81 +21,6 @@ module tally
contains
!===============================================================================
! CALCULATE_KEFF calculates the single cycle estimate of keff as well as the
! mean and standard deviation of the mean for active cycles and displays them
!===============================================================================
subroutine calculate_keff(i_cycle)
integer, intent(in) :: i_cycle ! index of current cycle
integer(8) :: total_bank ! total number of source sites
integer :: n ! active cycle number
real(8) :: k_cycle ! single cycle estimate of keff
real(8), save :: k_sum ! accumulated keff
real(8), save :: k_sum_sq ! accumulated keff**2
message = "Calculate cycle keff..."
call write_message(8)
! initialize sum and square of sum at beginning of run
if (i_cycle == 1) then
k_sum = ZERO
k_sum_sq = ZERO
end if
#ifdef MPI
! Collect number bank sites onto master process
call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, &
& MPI_COMM_WORLD, mpi_err)
#else
total_bank = n_bank
#endif
! Collect statistics and print output
if (master) then
! Since the creation of bank sites was originally weighted by the last
! cycle keff, we need to multiply by that keff to get the current cycle's
! value
k_cycle = real(total_bank)/real(n_particles)*keff
if (i_cycle > n_inactive) then
! Active cycle number
n = i_cycle - n_inactive
! Accumulate cycle estimate of k
k_sum = k_sum + k_cycle
k_sum_sq = k_sum_sq + k_cycle*k_cycle
! Determine mean and standard deviation of mean
keff = k_sum/n
keff_std = sqrt((k_sum_sq/n - keff*keff)/n)
! Display output for this cycle
if (i_cycle > n_inactive+1) then
write(UNIT=OUTPUT_UNIT, FMT=101) i_cycle, k_cycle, keff, keff_std
else
write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle
end if
else
! Display output for inactive cycle
write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle
keff = k_cycle
end if
end if
#ifdef MPI
! Broadcast new keff value to all processors
call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err)
#endif
100 format (2X,I5,2X,F8.5)
101 format (2X,I5,2X,F8.5,9X,F8.5,1X,F8.5)
end subroutine calculate_keff
!===============================================================================
! CREATE_TALLY_MAP creates a map that allows a quick determination of which
! tallies and bins need to be scored to when a particle makes a collision. This
@ -1043,7 +966,7 @@ contains
t => tallies(i)
! Write header block
call header("TALLY " // trim(int_to_str(t % id)), 3, UNIT_TALLY)
call header("TALLY " // trim(to_str(t % id)), unit=UNIT_TALLY, level=3)
! Handle surface current tallies separately
if (t % surface_current) then
@ -1125,8 +1048,8 @@ contains
do k = 1, t % n_macro_bins
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), macro_name(abs(t % macro_bins(k) % scalar)), &
real_to_str(t % scores(score_index,k) % val), &
trim(real_to_str(t % scores(score_index,k) % val_sq))
to_str(t % scores(score_index,k) % val), &
trim(to_str(t % scores(score_index,k) % val_sq))
end do
indent = indent - 2
@ -1176,14 +1099,14 @@ contains
end if
do i = 1, m % dimension(1)
string = "Mesh Index (" // trim(int_to_str(i)) // ", "
string = "Mesh Index (" // trim(to_str(i)) // ", "
len1 = len_trim(string)
do j = 1, m % dimension(2)
string = string(1:len1+1) // trim(int_to_str(j)) // ", "
string = string(1:len1+1) // trim(to_str(j)) // ", "
len2 = len_trim(string)
do k = 1, m % dimension(3)
! Write mesh cell index
string = string(1:len2+1) // trim(int_to_str(k)) // ")"
string = string(1:len2+1) // trim(to_str(k)) // ")"
write(UNIT=UNIT_TALLY, FMT='(1X,A)') trim(string)
do l = 1, n
@ -1202,15 +1125,15 @@ contains
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Left", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
bins(TS_SURFACE) = OUT_RIGHT
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Left", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
! Right Surface
bins(1:3) = (/ i, j, k /) + 1
@ -1218,15 +1141,15 @@ contains
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Right", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
bins(TS_SURFACE) = OUT_RIGHT
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Right", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
! Back Surface
bins(1:3) = (/ i, j-1, k /) + 1
@ -1234,15 +1157,15 @@ contains
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Back", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
bins(TS_SURFACE) = OUT_FRONT
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Back", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
! Front Surface
bins(1:3) = (/ i, j, k /) + 1
@ -1250,15 +1173,15 @@ contains
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Front", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
bins(TS_SURFACE) = OUT_FRONT
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Front", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
! Bottom Surface
bins(1:3) = (/ i, j, k-1 /) + 1
@ -1266,15 +1189,15 @@ contains
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Bottom", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
bins(TS_SURFACE) = OUT_TOP
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Bottom", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
! Top Surface
bins(1:3) = (/ i, j, k /) + 1
@ -1282,15 +1205,15 @@ contains
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Top", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
bins(TS_SURFACE) = OUT_TOP
score_index = sum((bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Top", &
real_to_str(t % scores(score_index,1) % val), &
trim(real_to_str(t % scores(score_index,1) % val_sq))
to_str(t % scores(score_index,1) % val), &
trim(to_str(t % scores(score_index,1) % val_sq))
end do
end do
@ -1320,38 +1243,38 @@ contains
select case(filter_type)
case (T_UNIVERSE)
i = t % universe_bins(bin) % scalar
label = int_to_str(universes(i) % id)
label = to_str(universes(i) % id)
case (T_MATERIAL)
i = t % material_bins(bin) % scalar
label = int_to_str(materials(i) % id)
label = to_str(materials(i) % id)
case (T_CELL)
i = t % cell_bins(bin) % scalar
label = int_to_str(cells(i) % id)
label = to_str(cells(i) % id)
case (T_CELLBORN)
i = t % cellborn_bins(bin) % scalar
label = int_to_str(cells(i) % id)
label = to_str(cells(i) % id)
case (T_SURFACE)
i = t % surface_bins(bin) % scalar
label = int_to_str(surfaces(i) % id)
label = to_str(surfaces(i) % id)
case (T_MESH)
m => meshes(t % mesh)
allocate(ijk(m % n_dimension))
call bin_to_mesh_indices(m, bin, ijk)
if (m % n_dimension == 2) then
label = "Index (" // trim(int_to_str(ijk(1))) // ", " // &
trim(int_to_str(ijk(2))) // ")"
label = "Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ")"
elseif (m % n_dimension == 3) then
label = "Index (" // trim(int_to_str(ijk(1))) // ", " // &
trim(int_to_str(ijk(2))) // ", " // trim(int_to_str(ijk(3))) // ")"
label = "Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
end if
case (T_ENERGYIN)
E0 = t % energy_in(bin)
E1 = t % energy_in(bin + 1)
label = "[" // trim(real_to_str(E0)) // ", " // trim(real_to_str(E1)) // ")"
label = "[" // trim(to_str(E0)) // ", " // trim(to_str(E1)) // ")"
case (T_ENERGYOUT)
E0 = t % energy_out(bin)
E1 = t % energy_out(bin + 1)
label = "[" // trim(real_to_str(E0)) // ", " // trim(real_to_str(E1)) // ")"
label = "[" // trim(to_str(E0)) // ", " // trim(to_str(E1)) // ")"
end select
end function get_label

4
src/xml-fortran/Makefile
View file

@ -8,6 +8,10 @@ objects = $(source:.f90=.o)
# Ignore unused variables
ifeq ($(MACHINE),bluegene)
override F90 = xlf2003
endif
ifeq ($(F90),ifort)
override F90FLAGS += -warn nounused
endif

1
src/xml-fortran/templates/settings_t.xml
View file

@ -26,5 +26,6 @@
<variable name="cutoff_" tag="cutoff" type="cutoff_xml" dimension="1" />
<variable name="cross_sections_" tag="cross_sections" type="word" length="255" />
<variable name="trace_" tag="trace" type="integer-array" />
<variable name="entropy_box_" tag="entropy_box" type="double-array" />
</template>