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

First version with tally capabilities -- still experimental.

Browse source
This commit is contained in:
Paul Romano 2011-04-26 03:16:42 +00:00
parent 620aed0013
commit 88dbd92f77
14 changed files with 734 additions and 138 deletions
Download patch file Download diff file Expand all files Collapse all files

21
ChangeLog
View file

@ -1,3 +1,24 @@
2011-04-25 Paul Romano <romano7@mit.edu>
* ace.f90: Fixed indentation of get_int and get_real. Added
subroutine get_macro_xs.
* fileio.f90: Added ability to read tally cards through read_tally
subroutine. Changed ERROR_CODE to ERROR_INT.
* global.f90: Added variables for tallies. Changed ERROR_CODE to
ERROR_INT and added ERROR_REAL.
* input_sample: Added tally card.
* main.f90: Use tallies_on variable when appropriate.
* output.f90: Add preliminary tally printing capabilites through
print_tally subroutine.
* physics.f90: Add call to score_tally for tallies. Changed
Sigma_t to Sigam_rxn.
* score.f90: Added score_tally and add_to_score subroutine.
* search.f90: Added interpolate subroutine.
* string.f90: Added is_number function.
* tallies.f90: CHANGED NAME TO score.f90
* types.f90: Added list of tallies to Universe, Lattice and
Cell. Modified Tally type.
2011-04-22 Paul Romano <romano7@mit.edu>
* geometry.f90: Deallocated count lists in neighbor_lists.

8
src/Makefile
View file

@ -17,7 +17,7 @@ modules = ace.f90 \
search.f90 \
source.f90 \
string.f90 \
tallies.f90 \
score.f90 \
types.f90
main_objects = $(modules:.f90=.o) $(main:.f90=.o)
@ -70,14 +70,14 @@ geometry.o: types.o global.o output.o string.o data_structures.o
global.o: types.o
main.o: global.o fileio.o output.o geometry.o mcnp_random.o \
source.o physics.o cross_section.o data_structures.o \
ace.o energy_grid.o tallies.o
ace.o energy_grid.o score.o
mpi_routines.o: global.o output.o types.o mcnp_random.o source.o
output.o: global.o types.o data_structures.o endf.o
physics.o: types.o global.o mcnp_random.o geometry.o output.o \
search.o endf.o
search.o endf.o score.o
score.o: global.o output.o ace.o
search.o: output.o
source.o: global.o mcnp_random.o output.o physics.o
string.o: global.o output.o
tallies.o: global.o output.o
types.o:
unittest.o: global.o energy_grid.o mpi_routines.o

93
src/ace.f90
View file

@ -1075,29 +1075,96 @@ contains
! in the XSS array
!=====================================================================
function get_int(n_values) result(array)
function get_int(n_values) result(array)
integer, intent(in) :: n_values ! number of values to read
integer :: array(n_values) ! array of values
integer, intent(in) :: n_values ! number of values to read
integer :: array(n_values) ! array of values
array = int(XSS(XSS_index:XSS_index + n_values - 1))
XSS_index = XSS_index + n_values
end function get_int
array = int(XSS(XSS_index:XSS_index + n_values - 1))
XSS_index = XSS_index + n_values
end function get_int
!=====================================================================
! GET_REAL returns an array of real(8)s read from the current position
! in the XSS array
!=====================================================================
function get_real(n_values) result(array)
function get_real(n_values) result(array)
integer, intent(in) :: n_values ! number of values to read
real(8) :: array(n_values) ! array of values
integer, intent(in) :: n_values ! number of values to read
real(8) :: array(n_values) ! array of values
array = XSS(XSS_index:XSS_index + n_values - 1)
XSS_index = XSS_index + n_values
array = XSS(XSS_index:XSS_index + n_values - 1)
XSS_index = XSS_index + n_values
end function get_real
end function get_real
!=====================================================================
! GET_MACRO_XS
!=====================================================================
function get_macro_xs(p, mat, MT) result(xs)
type(Particle), pointer :: p
type(Material), pointer :: mat
integer, intent(in) :: MT
real(8) :: xs
integer :: i
integer :: n_isotopes
integer :: IE
real(8) :: density
real(8) :: density_i
real(8) :: sigma_i
real(8) :: f
type(AceContinuous), pointer :: table => null()
type(AceReaction), pointer :: rxn => null()
! initialize xs
xs = ZERO
! find material atom density
density = cMaterial%atom_density
! loop over all isotopes in material
n_isotopes = cMaterial % n_isotopes
do i = 1, n_isotopes
table => xs_continuous(cMaterial % table(i))
! determine nuclide atom density
density_i = cMaterial%atom_percent(i) * density
! search nuclide energy grid
IE = table%grid_index(p % IE)
f = (p%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE))
! handle special case of total cross section
if (MT == 1) then
xs = xs + density * (ONE-f) * table%sigma_t(IE) + &
& f * (table%sigma_t(IE+1))
cycle
end if
! loop over reactions in isotope
do j = 1, table % n_reaction
rxn => table % reactions(i)
! check for matching MT
if (MT /= rxn % MT) cycle
! if energy is below threshold for this reaction, skip it
if (IE < rxn % IE) cycle
! add to cumulative probability
sigma_i = (ONE-f) * rxn%sigma(IE-rxn%IE+1) + &
& f * (rxn%sigma(IE-rxn%IE+2))
end do
! calculate nuclide macroscopic cross-section
xs = xs + density_i * sigma_i
end do
end function get_macro_xs
end module ace

197
src/fileio.f90
View file

@ -2,7 +2,8 @@ module fileio
use global
use types, only: Cell, Surface, ExtSource, ListKeyValueCI
use string, only: split_string_wl, lower_case, split_string, concatenate
use string, only: split_string_wl, lower_case, split_string, &
& concatenate, is_number
use output, only: message, warning, error
use data_structures, only: dict_create, dict_add_key, dict_get_key, &
& dict_has_key, DICT_NULL, dict_keys, list_size
@ -12,7 +13,6 @@ module fileio
type(DictionaryII), pointer :: & ! used to count how many cells each
& ucount_dict => null(), & ! universe contains
& bc_dict => null() ! store boundary conditions
integer, allocatable :: index_cell_in_univ(:)
@ -93,6 +93,7 @@ contains
n_universes = 0
n_surfaces = 0
n_materials = 0
n_tallies = 0
call dict_create(ucount_dict)
call dict_create(universe_dict)
@ -131,6 +132,8 @@ contains
n_surfaces = n_surfaces + 1
case ('material')
n_materials = n_materials + 1
case ('tally')
n_tallies = n_tallies + 1
case ('lattice')
n_lattices = n_lattices + 1
@ -167,6 +170,7 @@ contains
allocate(lattices(n_lattices))
allocate(surfaces(n_surfaces))
allocate(materials(n_materials))
allocate(tallies(n_tallies))
! Also allocate a list for keeping track of where cells have been
! assigned in each universe
@ -178,6 +182,7 @@ contains
call dict_create(surface_dict)
call dict_create(bc_dict)
call dict_create(material_dict)
call dict_create(tally_dict)
! We also need to allocate the cell count lists for each
! universe. The logic for this is a little more convoluted. In
@ -227,6 +232,7 @@ contains
integer :: index_lattice ! index in lattices array
integer :: index_material ! index in materials array
integer :: index_source ! index in source array (?)
integer :: index_tally ! index in tally array
character(250) :: line ! a line of words
character(250) :: msg ! output/error message
character(32) :: words(max_words) ! words on a single line
@ -239,6 +245,7 @@ contains
index_surface = 0
index_lattice = 0
index_material = 0
index_tally = 0
index_source = 0
open(FILE=filename, UNIT=in, STATUS='old', &
@ -280,6 +287,11 @@ contains
index_material = index_material + 1
call read_material(index_material, words, n)
case ('tally')
! Read tally entry
index_tally = index_tally + 1
call read_tally(index_tally, words, n)
case ('source')
call read_source(words, n)
@ -291,7 +303,7 @@ contains
case ('verbosity')
verbosity = str_to_int(words(2))
if (verbosity == ERROR_CODE) then
if (verbosity == ERROR_INT) then
msg = "Invalid verbosity: " // words(2)
call error(msg)
end if
@ -468,7 +480,7 @@ contains
! Read cell identifier
c % uid = str_to_int(words(2))
if (c % uid == ERROR_CODE) then
if (c % uid == ERROR_INT) then
msg = "Invalid cell name: " // words(2)
call error(msg)
end if
@ -476,7 +488,7 @@ contains
! Read cell universe
universe_num = str_to_int(words(3))
if (universe_num == ERROR_CODE) then
if (universe_num == ERROR_INT) then
msg = "Invalid universe: " // words(3)
call error(msg)
end if
@ -489,7 +501,7 @@ contains
! find universe
universe_num = str_to_int(words(5))
if (universe_num == ERROR_CODE) then
if (universe_num == ERROR_INT) then
msg = "Invalid universe fill: " // words(5)
call error(msg)
end if
@ -506,7 +518,7 @@ contains
c % type = CELL_NORMAL
c % material = str_to_int(words(4))
c % fill = 0
if (c % material == ERROR_CODE) then
if (c % material == ERROR_INT) then
msg = "Invalid material number: " // words(4)
call error(msg)
end if
@ -647,7 +659,7 @@ contains
! Read surface identifier
surface_uid = str_to_int(words(2))
if (surface_uid == ERROR_CODE) then
if (surface_uid == ERROR_INT) then
msg = "Invalid surface name: " // words(2)
call error(msg)
end if
@ -695,7 +707,7 @@ contains
! Read lattice universe
universe_num = str_to_int(words(2))
if (universe_num == ERROR_CODE) then
if (universe_num == ERROR_INT) then
msg = "Invalid universe: " // words(2)
call error(msg)
end if
@ -717,10 +729,10 @@ contains
! Read number of lattice cells in each direction
n_x = str_to_int(words(4))
n_y = str_to_int(words(5))
if (n_x == ERROR_CODE) then
if (n_x == ERROR_INT) then
msg = "Invalid number of lattice cells in x-direction: " // words(4)
call error(msg)
elseif (n_y == ERROR_CODE) then
elseif (n_y == ERROR_INT) then
msg = "Invalid number of lattice cells in y-direction: " // words(5)
call error(msg)
end if
@ -748,7 +760,7 @@ contains
do i = 1, n_x
index_word = 9 + j*n_x + i
universe_num = str_to_int(words(index_word))
if (universe_num == ERROR_CODE) then
if (universe_num == ERROR_INT) then
msg = "Invalid universe number: " // words(index_word)
call error(msg)
end if
@ -799,6 +811,161 @@ contains
end subroutine read_source
!=====================================================================
! READ_TALLY
!=====================================================================
subroutine read_tally(index, words, n_words)
integer, intent(in) :: index ! index in materials array
character(*), intent(in) :: words(n_words) ! words on material entry
integer, intent(in) :: n_words ! number of words
integer :: i ! index in words array
integer :: j
integer :: k
integer :: count
integer :: MT
integer :: cell_uid
integer :: r_bins, c_bins, e_bins
real(8) :: E
character(32) :: word
character(250) :: msg
type(Tally), pointer :: t => null()
t => tallies(index)
! Read tally identifier
t % uid = str_to_int(words(2))
if (t % uid == ERROR_INT) then
msg = "Invalid tally name: " // words(2)
call error(msg)
end if
call dict_add_key(tally_dict, t % uid, index)
i = 3
do while (i <= n_words)
word = words(i)
call lower_case(word)
select case (trim(word))
case ('reaction')
! ==========================================================
! READ REACTION LIST
! Determine how many reactions are listed
do j = i+1, n_words
if (.not. is_number(words(j))) exit
end do
count = j - (i+1)
! Allocate space to store reactions
allocate(t % reactions(count))
! Read reaction MT values
do j = 1, count
MT = str_to_int(words(i+j))
if (MT == ERROR_INT) then
msg = "Invalid reaction MT on tally: " // int_to_str(t%uid)
call error(msg)
end if
t % reactions(j) = MT
end do
! Determine whether to sum reactions or individually bin
if (trim(words(i+count+1)) == 'sum') then
t % reaction_type = TALLY_SUM
r_bins = 1
else
t % reaction_type = TALLY_BINS
r_bins = count
end if
case ('cell')
! ==========================================================
! READ CELL LIST
! Determine how many reactions are listed
do j = i+1, n_words
if (.not. is_number(words(j))) exit
end do
count = j - (i+1)
! Allocate space to store reactions
allocate(t % cells(count))
! Read reaction MT values
do j = 1, count
cell_uid = str_to_int(words(i+j))
if (cell_uid == ERROR_INT) then
msg = "Invalid cell number on tally: " // int_to_str(t%uid)
call error(msg)
end if
t % cells(j) = cell_uid
end do
! Determine whether to sum reactions or individually bin
if (trim(words(i+count+1)) == 'sum') then
t % cell_type = TALLY_SUM
c_bins = 1
else
t % cell_type = TALLY_BINS
c_bins = count
end if
case ('energy')
! ==========================================================
! READ ENERGY LIST
! Determine how many energies are listed
do j = i+1, n_words
if (.not. is_number(words(j)(1:1))) exit
end do
count = j - (i+1)
! Allocate space to store energies
allocate(t % energies(count))
! Read reaction MT values
do j = 1, count
E = str_to_real(words(i+j))
if (E == ERROR_REAL) then
msg = "Invalid energy on tally: " // int_to_str(t % uid)
call error(msg)
end if
t % energies(j) = E
end do
e_bins = count - 1
case ('material')
! TODO: Add ability to read material lists
case ('universe')
! TODO: Add ability to read universe lists
case ('lattice')
! TODO: Add lattice tally ability
end select
! Move index in words forward
i = i + 1 + count
end do
! allocate tally scores
allocate(t % score(r_bins, c_bins, e_bins))
! initialize tallies
forall (i=1:r_bins, j=1:c_bins, k=1:e_bins)
t % score(i,j,k) % n_events = 0
t % score(i,j,k) % val = ZERO
t % score(i,j,k) % val_sq = ZERO
end forall
end subroutine read_tally
!=====================================================================
! READ_MATERIAL parses a material card. Note that atom percents and
! densities are normalized in a separate routine
@ -964,21 +1131,21 @@ contains
! Read number of cycles
n_cycles = str_to_int(words(2))
if (n_cycles == ERROR_CODE) then
if (n_cycles == ERROR_INT) then
msg = "Invalid number of cycles: " // words(2)
call error(msg)
end if
! Read number of inactive cycles
n_inactive = str_to_int(words(3))
if (n_inactive == ERROR_CODE) then
if (n_inactive == ERROR_INT) then
msg = "Invalid number of inactive cycles: " // words(2)
call error(msg)
end if
! Read number of particles
n_particles = str_to_int(words(4))
if (n_particles == ERROR_CODE) then
if (n_particles == ERROR_INT) then
msg = "Invalid number of particles: " // words(2)
call error(msg)
end if

35
src/global.f90
View file

@ -16,11 +16,15 @@ module global
type(Material), allocatable, target :: materials(:)
type(Isotope), allocatable, target :: isotopes(:)
type(xsData), allocatable, target :: xsdatas(:)
integer :: n_cells ! # of cells
integer :: n_universes ! # of universes
integer :: n_lattices ! # of lattices
integer :: n_surfaces ! # of surfaces
integer :: n_materials ! # of materials
type(Tally), allocatable, target :: tallies(:)
type(Tally), allocatable, target :: tallies_global(:)
integer :: n_cells ! # of cells
integer :: n_universes ! # of universes
integer :: n_lattices ! # of lattices
integer :: n_surfaces ! # of surfaces
integer :: n_materials ! # of materials
integer :: n_tallies ! # of tallies
integer :: n_tallies_global ! # of global tallies
! These dictionaries provide a fast lookup mechanism
type(DictionaryII), pointer :: cell_dict
@ -29,6 +33,7 @@ module global
type(DictionaryII), pointer :: surface_dict
type(DictionaryII), pointer :: material_dict
type(DictionaryII), pointer :: isotope_dict
type(DictionaryII), pointer :: tally_dict
type(DictionaryCI), pointer :: xsdata_dict
type(DictionaryCI), pointer :: ace_dict
@ -68,6 +73,8 @@ module global
! cycle keff
real(8) :: keff
logical :: tallies_on
! Parallel processing variables
integer :: n_procs ! number of processes
integer :: rank ! rank of process
@ -223,6 +230,13 @@ module global
& N_PT = 116, &
& N_DA = 117
! Tally distinguishers
integer, parameter :: &
& TALLY_FLUX = -1, & ! tally the flux only
& TALLY_ALL = 0, & ! tally all reactions
& TALLY_BINS = 1, & ! tally individual (reactions, cells)
& TALLY_SUM = 2 ! tally sum of specified (reactions, cells)
! Fission neutron emission (nu) type
integer, parameter :: &
& NU_NONE = 0, & ! No nu values (non-fissionable)
@ -232,7 +246,10 @@ module global
! Integer code for read error -- better hope this number is never
! used in an input file!
integer, parameter :: &
& ERROR_CODE = -huge(0)
& ERROR_INT = -huge(0)
real(8), parameter :: &
& ERROR_REAL = -huge(0.0_8) * 0.917826354_8
! The verbosity controls how much information will be printed to the
! screen and in logs
@ -358,7 +375,7 @@ contains
! read string into integer
read(UNIT=str, FMT=fmt, IOSTAT=ioError) num
if (ioError > 0) num = ERROR_CODE
if (ioError > 0) num = ERROR_INT
end function str_to_int
@ -409,9 +426,7 @@ contains
! Read string
read(string, fmt, iostat=readError) num
if (readError > 0) then
! return error code for real?
end if
if (readError > 0) num = ERROR_REAL
end function str_to_real

2
src/input_sample
View file

@ -20,6 +20,8 @@ material 41 -1.0 &
source box -4 -4 -4 4 4 4
tally 1 reaction 2 cell 2 energy limits
xs_data /opt/serpent/xsdata/endfb7
criticality 15 5 10000
verbosity 7

9
src/main.f90
View file

@ -13,7 +13,7 @@ program main
use ace, only: read_xs
use energy_grid, only: unionized_grid, original_indices
use mpi_routines, only: setup_mpi, synchronize_bank, t_sync
use tallies, only: calculate_keff
use score, only: calculate_keff
#ifdef MPI
use mpi
@ -111,6 +111,8 @@ contains
msg = "Running problem..."
call message(msg, 6)
tallies_on = .false.
#ifdef MPI
t0 = MPI_WTIME()
#endif
@ -148,9 +150,10 @@ contains
! print cycle information
end do CYCLE_LOOP
! Turn tallies on once inactive cycles are complete
if (i_cycle == n_inactive) tallies_on = .true.
! Collect all tallies and print
end do CYCLE_LOOP
#ifdef MPI
! print run time

74
src/output.f90
View file

@ -511,6 +511,66 @@ contains
end subroutine print_material
!=====================================================================
! PRINT_TALLY displays the attributes of a tally
!=====================================================================
subroutine print_tally(tal)
type(Tally), pointer :: tal
integer :: i
integer :: MT
character(250) :: string
write(ou,*) 'Tally ' // int_to_str(tal % uid)
select case (tal % reaction_type)
case (TALLY_FLUX)
write(ou,*) ' Type: Flux'
case (TALLY_ALL)
write(ou,*) ' Type: Total Collision Rate'
case (TALLY_BINS)
write(ou,*) ' Type: Partial reactions'
case (TALLY_SUM)
write(ou,*) ' Type: Partial reactions (summed)'
end select
select case (tal % cell_type)
case (TALLY_BINS)
write(ou,*) ' Cell Type: Separate bins'
case (TALLY_SUM)
write(ou,*) ' Cell Type: Sum over cells'
end select
if (allocated(tal % reactions)) then
string = ""
do i = 1, size(tal % reactions)
MT = tal % reactions(i)
string = trim(string) // ' ' // trim(reaction_name(MT))
end do
write(ou,*) ' Reactions:' // trim(string)
end if
if (allocated(cells)) then
string = ""
do i = 1, size(tal % cells)
string = trim(string) // ' ' // trim(int_to_str(tal % cells(i)))
end do
write(ou,*) ' Cells:' // trim(string)
end if
if (allocated(tal % energies)) then
string = ""
do i = 1, size(tal % energies)
string = trim(string) // ' ' // trim(real_to_str(tal % energies(i)))
end do
write(ou,*) ' Energies:' // trim(string)
end if
write(ou,*)
end subroutine print_tally
!=====================================================================
! PRINT_SUMMARY displays the attributes of all cells, universes,
! surfaces and materials read in the input file. Very useful for
@ -524,6 +584,7 @@ contains
type(Universe), pointer :: u => null()
type(Lattice), pointer :: l => null()
type(Material), pointer :: m => null()
type(Tally), pointer :: t => null()
integer :: i
! print summary of cells
@ -578,11 +639,24 @@ contains
call print_material(m)
end do
! print summary of tallies
if (n_tallies > 0) then
write(ou,*) '============================================='
write(ou,*) '=> TALLY SUMMARY <='
write(ou,*) '============================================='
write(ou,*)
do i = 1, n_tallies
t=> tallies(i)
call print_tally(t)
end do
end if
nullify(s)
nullify(c)
nullify(u)
nullify(l)
nullify(m)
nullify(t)
end subroutine print_summary

32
src/physics.f90
View file

@ -8,6 +8,7 @@ module physics
use output, only: error, warning, message, print_particle
use search, only: binary_search
use endf, only: reaction_name
use score, only: score_tally
implicit none
@ -81,8 +82,6 @@ contains
p%xyz = p%xyz + distance * p%uvw
p%xyz_local = p%xyz_local + distance * p%uvw
! Add pathlength tallies
if (d_to_collision > d_to_boundary) then
p % cell = 0
if (in_lattice) then
@ -159,50 +158,51 @@ contains
real(8) :: density_i ! atom density of nuclide
real(8) :: prob ! cumulative probability
real(8) :: r1 ! random number
real(8), allocatable :: Sigma_t(:) ! macroscopic xs for each nuclide
real(8) :: flux ! collision estimator of flux
real(8), allocatable :: Sigma_rxn(:) ! macroscopic xs for each nuclide
character(250) :: msg ! output/error message
type(AceContinuous), pointer :: table
type(AceReaction), pointer :: rxn
! TODO: tallies
density = cMaterial%atom_density
! calculate total cross-section for each nuclide at current energy
! in order to create discrete pdf for sampling nuclide
n_isotopes = cMaterial%n_isotopes
allocate(Sigma_t(n_isotopes))
allocate(Sigma_rxn(n_isotopes))
do i = 1, n_isotopes
table => xs_continuous(cMaterial%table(i))
density_i = cMaterial%atom_percent(i)*density
! determine nuclide atom density
density_i = cMaterial%atom_percent(i) * density
! search nuclide energy grid
IE = table%grid_index(p % IE)
f = (p%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE))
! calculate macroscopic xs for this nuclide
Sigma_t(i) = density_i*((ONE-f)*table%Sigma_t(IE) + &
! calculate nuclide macroscopic cross-section
Sigma_rxn(i) = density_i*((ONE-f)*table%Sigma_t(IE) + &
& f*(table%Sigma_t(IE+1)))
end do
total = sum(Sigma_t)
total = sum(Sigma_rxn)
! sample nuclide
r1 = rang()
prob = ZERO
do i = 1, n_isotopes
prob = prob + Sigma_t(i) / total
prob = prob + Sigma_rxn(i) / total
if (r1 < prob) exit
end do
! Get table, total xs, interpolation factor
density_i = cMaterial%atom_percent(i)*density
table => xs_continuous(cMaterial%table(i))
sigma = Sigma_t(i) / density_i
sigma = Sigma_rxn(i) / density_i
IE = table%grid_index(p % IE)
f = (p%E - table%energy(IE))/(table%energy(IE+1) - table%energy(IE))
! free memory
deallocate(Sigma_t)
deallocate(Sigma_rxn)
! sample reaction channel
r1 = rang()*sigma
@ -248,6 +248,12 @@ contains
call warning(msg)
end select
! Add collision estimator tallies
if (tallies_on) then
flux = p%wgt / SIGMA
call score_tally(p, flux)
end if
! check for very low energy
if (p % E < 1.0e-100_8) then
p % alive = .false.

277
src/score.f90 Normal file
View file

@ -0,0 +1,277 @@
module score
use global
use output, only: message, error
use ace, only: get_macro_xs
use search, only: binary_search
#ifdef MPI
use mpi
#endif
implicit none
contains
!=====================================================================
! CALCULATE_KEFF
!=====================================================================
subroutine calculate_keff(i_cycle)
integer, intent(in) :: i_cycle ! index of current cycle
integer(8) :: total_bank
integer :: n
integer :: ierr
real(8) :: kcoll ! keff collision estimator
real(8) :: ktemp ! MPI-reduced keff and stdev
real(8), save :: k1 = 0. ! accumulated keff
real(8), save :: k2 = 0. ! accumulated keff**2
real(8) :: std ! stdev of keff over active cycles
character(250) :: msg
msg = "Calculate cycle keff..."
call message(msg, 8)
! set k1 and k2 at beginning of run
if (i_cycle == 1) then
k1 = ZERO
k2 = 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, ierr)
#else
total_bank = n_bank
#endif
! Collect statistics and print output
if (master) then
kcoll = real(total_bank)/real(n_particles)*keff
if (i_cycle > n_inactive) then
n = i_cycle - n_inactive
k1 = k1 + kcoll
k2 = k2 + kcoll**2
keff = k1/n
std = sqrt((k2/n-keff**2)/n)
if (i_cycle > n_inactive+1) then
write(6,101) i_cycle, kcoll, keff, std
else
write(6,100) i_cycle, kcoll
end if
else
write(6,100) i_cycle, kcoll
keff = kcoll
end if
end if
#ifdef MPI
call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, ierr)
#endif
100 format (2X,I4,3X,F7.5)
101 format (2X,I4,3X,F7.5,10X,F7.5,2X,F7.5)
end subroutine calculate_keff
!=====================================================================
! SCORE_TALLY
!=====================================================================
subroutine score_tally(p, flux)
type(Particle), pointer :: p ! particle
real(8), intent(in) :: flux ! estimator of flux
integer :: i ! index for tallies in cell
integer :: j ! index for reactions in tally
integer :: e_bin ! energy bin
integer :: c_bin ! cell bin
integer :: r_bin ! reaction bin
integer :: n_energy ! number of energy bins
integer :: n_reaction ! number of reactions
integer :: type ! cell or reaction type
integer :: MT ! reaction MT value
real(8) :: E ! energy of particle
real(8) :: val ! value to score
real(8) :: Sigma ! macroscopic cross section of reaction
character(250) :: msg ! output/error message
type(Cell), pointer :: c => null()
type(Tally), pointer :: t => null()
! ================================================================
! HANDLE LOCAL TALLIES
c => cells(p % cell)
if (.not. allocated(c % tallies)) return
! if so, loop over each tally
do i = 1, size(c % tallies)
t => tallies(c % tallies(i))
! =============================================================
! DETERMINE ENERGY BIN
if (allocated(t % energies)) then
E = p % E
n_energy = size(t % energies)
if (E < t % energies(1) .or. E > t % energies(n_energy)) then
! Energy outside of specified range
cycle
else
e_bin = binary_search(t % energies, n_energy, E)
end if
end if
! =============================================================
! DETERMINE CELL BIN
type = t % cell_type
if (type == TALLY_SUM) then
! Sum tallies from separate cells into one bin
c_bin = 1
elseif (type == TALLY_BINS) then
! Need to determine cell bin
do j = 1, size(t % cells)
if (p % cell == t % cells(j)) then
c_bin = j
exit
end if
end do
else
msg = "Invalid type for cell bins in tally " // int_to_str(t % uid)
call error(msg)
end if
! =============================================================
! DETERMINE REACTION BIN AND ADD VALUE TO SCORE
type = t % reaction_type
if (type == TALLY_FLUX) then
! Tally flux only
val = flux
call add_to_score(t % score(r_bin, c_bin, e_bin), val)
elseif (type == TALLY_ALL) then
! Tally total reaction rate
val = ONE
r_bin = 1
call add_to_score(t % score(r_bin, c_bin, e_bin), val)
elseif (type == TALLY_BINS) then
! Individually bin reactions
n_reaction = t % reaction_type
r_bin = 0
do j = 1, n_reaction
MT = t % reactions(j)
Sigma = get_macro_xs(p, cMaterial, MT)
val = Sigma * flux
r_bin = r_bin + 1
call add_to_score(t % score(r_bin, c_bin, e_bin), &
& val)
end do
elseif (type == TALLY_SUM) then
! Tally reactions in one bin
n_reaction = t % reaction_type
r_bin = 1
do j = 1, n_reaction
MT = t % reactions(j)
Sigma = get_macro_xs(p, cMaterial, MT)
val = Sigma * flux
call add_to_score(t % score(r_bin, c_bin, e_bin), &
& val)
end do
end if
end do
! ================================================================
! HANDLE GLOBAL TALLIES
do i = 1, n_tallies_global
t => tallies_global(i)
! =============================================================
! DETERMINE ENERGY BIN
if (allocated(t % energies)) then
E = p % E
n_energy = size(t % energies)
if (E < t % energies(1) .or. E > t % energies(n_energy)) then
! Energy outside of specified range
cycle
else
e_bin = binary_search(t % energies, n_energy, E)
end if
end if
! Since it's a global tally, the cell bin is unity
c_bin = 1
! =============================================================
! DETERMINE REACTION BIN AND ADD VALUE TO SCORE
type = t % reaction_type
if (type == TALLY_FLUX) then
! Tally flux only
val = flux
call add_to_score(t % score(r_bin, c_bin, e_bin), val)
elseif (type == TALLY_ALL) then
! Tally total reaction rate
val = ONE
r_bin = 1
call add_to_score(t % score(r_bin, c_bin, e_bin), val)
elseif (type == TALLY_BINS) then
! Individually bin reactions
n_reaction = t % reaction_type
r_bin = 0
do j = 1, n_reaction
MT = t % reactions(j)
Sigma = get_macro_xs(p, cMaterial, MT)
val = Sigma * flux
r_bin = r_bin + 1
call add_to_score(t % score(r_bin, c_bin, e_bin), &
& val)
end do
elseif (type == TALLY_SUM) then
! Tally reactions in one bin
n_reaction = t % reaction_type
r_bin = 1
do j = 1, n_reaction
MT = t % reactions(j)
Sigma = get_macro_xs(p, cMaterial, MT)
val = Sigma * flux
call add_to_score(t % score(r_bin, c_bin, e_bin), &
& val)
end do
end if
end do
! ================================================================
! TODO: Add lattice tallies
! ================================================================
! TODO: Add mesh tallies
end subroutine score_tally
!=====================================================================
! ADD_TO_SCORE
!=====================================================================
subroutine add_to_score(score, val)
type(TallyScore), intent(inout) :: score
real(8), intent(in) :: val
score % n_events = score % n_events + 1
score % val = score % val + val
score % val_sq = score % val_sq + val*val
end subroutine add_to_score
end module score

16
src/search.f90
View file

@ -55,4 +55,20 @@ contains
end function binary_search
!=====================================================================
! INTERPOLATE
!=====================================================================
function interpolate(array, n, index, f) result(val)
real(8), intent(in) :: array(n)
integer, intent(in) :: n
integer, intent(in) :: index
real(8), intent(in) :: f
real(8) :: val
val = (ONE-f) * array(index) + f * array(index+1)
end function interpolate
end module search

22
src/string.f90
View file

@ -150,7 +150,6 @@ contains
!=====================================================================
elemental subroutine lower_case(word)
! convert a word to lower case
character(*), intent(inout) :: word
@ -164,4 +163,25 @@ contains
end subroutine lower_case
!=====================================================================
! IS_NUMBER determines whether a string of characters is all 0-9
! characters
!=====================================================================
function is_number(word) result(number)
character(*), intent(in) :: word
logical :: number
integer :: i
integer :: ic
number = .true.
do i = 1, len_trim(word)
ic = ichar(word(i:i))
if (ic < 48 .or. ic >= 58) number = .false.
end do
end function is_number
end module string

78
src/tallies.f90
View file

@ -1,78 +0,0 @@
module tallies
use global
use output, only: message, error
#ifdef MPI
use mpi
#endif
implicit none
contains
!=====================================================================
! CALCULATE_KEFF
!=====================================================================
subroutine calculate_keff(i_cycle)
integer, intent(in) :: i_cycle ! index of current cycle
integer(8) :: total_bank
integer :: n
integer :: ierr
real(8) :: kcoll ! keff collision estimator
real(8) :: ktemp ! MPI-reduced keff and stdev
real(8), save :: k1 = 0. ! accumulated keff
real(8), save :: k2 = 0. ! accumulated keff**2
real(8) :: std ! stdev of keff over active cycles
character(250) :: msg
msg = "Calculate cycle keff..."
call message(msg, 8)
! set k1 and k2 at beginning of run
if (i_cycle == 1) then
k1 = ZERO
k2 = 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, ierr)
#else
total_bank = n_bank
#endif
! Collect statistics and print output
if (master) then
kcoll = real(total_bank)/real(n_particles)*keff
if (i_cycle > n_inactive) then
n = i_cycle - n_inactive
k1 = k1 + kcoll
k2 = k2 + kcoll**2
keff = k1/n
std = sqrt((k2/n-keff**2)/n)
if (i_cycle > n_inactive+1) then
write(6,101) i_cycle, kcoll, keff, std
else
write(6,100) i_cycle, kcoll
end if
else
write(6,100) i_cycle, kcoll
keff = kcoll
end if
end if
#ifdef MPI
call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, ierr)
#endif
100 format (2X,I4,3X,F7.5)
101 format (2X,I4,3X,F7.5,10X,F7.5,2X,F7.5)
end subroutine calculate_keff
end module tallies

8
src/types.f90
View file

@ -15,6 +15,7 @@ module types
real(8) :: x0 ! Translation in x-coordinate
real(8) :: y0 ! Translation in y-coordinate
real(8) :: z0 ! Translation in z-coordinate
integer, allocatable :: tallies(:)
end type Universe
!=====================================================================
@ -33,6 +34,7 @@ module types
real(8) :: width_x ! width of lattice cell
real(8) :: width_y ! width of lattice cell
integer, allocatable :: element(:,:) ! specified universes
integer, allocatable :: tallies(:)
end type Lattice
!=====================================================================
@ -66,6 +68,7 @@ module types
& surfaces(:) ! List of surfaces bounding cell -- note
! that parentheses, union, etc operators
! will be listed here too
integer, allocatable :: tallies(:)
end type Cell
!=====================================================================
@ -124,6 +127,9 @@ module types
integer :: uid
integer :: type
real(8) :: volume
integer :: cell_type
integer :: reaction_type
integer :: material_type
integer, allocatable :: reactions(:)
integer, allocatable :: cells(:)
integer, allocatable :: materials(:)
@ -134,7 +140,7 @@ module types
integer :: n_x
integer :: n_y
integer :: n_z
type(TallyScore), allocatable :: score(:)
type(TallyScore), allocatable :: score(:,:,:)
end type Tally
!=====================================================================