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

Merge branches 'master' and 'cmfd' into cmfd

Browse source
This commit is contained in:
Bryan Herman 2012-01-28 15:57:58 -05:00
commit 254aa57ef9
11 changed files with 1018 additions and 530 deletions
Download patch file Download diff file Expand all files Collapse all files

80
src/constants.F90
View file

@ -224,47 +224,65 @@ module constants
! ============================================================================
! TALLY-RELATED CONSTANTS
! Tally macro reactions
integer, parameter :: N_MACRO_TYPES = 15
! Tally type
integer, parameter :: &
MACRO_FLUX = -1, & ! flux
MACRO_TOTAL = -2, & ! total reaction rate
MACRO_SCATTER = -3, & ! scattering rate
MACRO_NU_SCATTER = -4, & ! scattering production rate
MACRO_SCATTER_1 = -5, & ! first scattering moment
MACRO_SCATTER_2 = -6, & ! second scattering moment
MACRO_SCATTER_3 = -7, & ! third scattering moment
MACRO_N_1N = -8, & ! (n,1n) rate
MACRO_N_2N = -9, & ! (n,2n) rate
MACRO_N_3N = -10, & ! (n,3n) rate
MACRO_N_4N = -11, & ! (n,4n) rate
MACRO_ABSORPTION = -12, & ! absorption rate
MACRO_FISSION = -13, & ! fission rate
MACRO_NU_FISSION = -14, & ! neutron production rate
MACRO_CURRENT = -15 ! partial current
TALLY_VOLUME = 1, &
TALLY_SURFACE_CURRENT = 2
! Tally estimator types
integer, parameter :: &
ESTIMATOR_ANALOG = 1, &
ESTIMATOR_TRACKLENGTH = 2
! Event types for tallies
integer, parameter :: &
EVENT_SURFACE = -2, &
EVENT_LATTICE = -1, &
EVENT_SCATTER = 1, &
EVENT_ABSORB = 2, &
EVENT_FISSION = 3
! Tally score type
integer, parameter :: N_SCORE_TYPES = 15
integer, parameter :: &
SCORE_FLUX = -1, & ! flux
SCORE_TOTAL = -2, & ! total reaction rate
SCORE_SCATTER = -3, & ! scattering rate
SCORE_NU_SCATTER = -4, & ! scattering production rate
SCORE_SCATTER_1 = -5, & ! first scattering moment
SCORE_SCATTER_2 = -6, & ! second scattering moment
SCORE_SCATTER_3 = -7, & ! third scattering moment
SCORE_N_1N = -8, & ! (n,1n) rate
SCORE_N_2N = -9, & ! (n,2n) rate
SCORE_N_3N = -10, & ! (n,3n) rate
SCORE_N_4N = -11, & ! (n,4n) rate
SCORE_ABSORPTION = -12, & ! absorption rate
SCORE_FISSION = -13, & ! fission rate
SCORE_NU_FISSION = -14, & ! neutron production rate
SCORE_CURRENT = -15 ! partial current
! Tally map bin finding
integer, parameter :: NO_BIN_FOUND = -1
! Tally filter and map types
integer, parameter :: TALLY_TYPES = 8
integer, parameter :: N_FILTER_TYPES = 8
integer, parameter :: &
T_UNIVERSE = 1, &
T_MATERIAL = 2, &
T_CELL = 3, &
T_CELLBORN = 4, &
T_SURFACE = 5, &
T_MESH = 6, &
T_ENERGYIN = 7, &
T_ENERGYOUT = 8
FILTER_UNIVERSE = 1, &
FILTER_MATERIAL = 2, &
FILTER_CELL = 3, &
FILTER_CELLBORN = 4, &
FILTER_SURFACE = 5, &
FILTER_MESH = 6, &
FILTER_ENERGYIN = 7, &
FILTER_ENERGYOUT = 8
! Filter types for surface current tallies
integer, parameter :: &
TS_MESH_X = 1, &
TS_MESH_Y = 2, &
TS_MESH_Z = 3, &
TS_ENERGYIN = 4, &
TS_SURFACE = 5
SURF_FILTER_MESH_X = 1, &
SURF_FILTER_MESH_Y = 2, &
SURF_FILTER_MESH_Z = 3, &
SURF_FILTER_ENERGYIN = 4, &
SURF_FILTER_SURFACE = 5
! Tally surface current directions
integer, parameter :: &

194
src/hdf5_interface.F90
View file

@ -110,6 +110,9 @@ contains
call hdf5_write_geometry()
call hdf5_write_materials()
if (n_tallies > 0) then
call hdf5_write_tallies()
end if
end subroutine hdf5_write_summary
@ -422,6 +425,197 @@ contains
end subroutine hdf5_write_materials
!===============================================================================
! HDF5_WRITE_TALLIES
!===============================================================================
subroutine hdf5_write_tallies()
integer :: i, j
integer(SIZE_T) :: num_elements
integer(HSIZE_T) :: dims(1)
integer(HSIZE_T) :: coord(1,1)
integer(HSIZE_T) :: size_string = 12
integer(HID_T) :: string_type
integer(HID_T) :: tallies_group
integer(HID_T) :: temp_group
integer(HID_T) :: dspace
integer(HID_T) :: subspace
integer(HID_T) :: dset
character(12) :: string
type(TallyObject), pointer :: t => null()
type(StructuredMesh), pointer :: m => null()
! Create group for tallies
call h5gcreate_f(hdf5_output_file, "/tallies", tallies_group, hdf5_err)
! Use H5LT interface to write number of tallies
call hdf5_make_integer(tallies_group, "n_tallies", n_tallies)
! Write information on each material
do i = 1, n_tallies
t => tallies(i)
! Create group for i-th universe
call h5gcreate_f(tallies_group, "tally " // trim(to_str(t % id)), &
temp_group, hdf5_err)
! =======================================================================
! WRITE INFORMATION ON TALLY FILTERS
! Write filters
call hdf5_make_integer(temp_group, "n_filters", t % n_filters)
if (t % n_filters > 0) then
dims(1) = t % n_filters
call h5ltmake_dataset_int_f(temp_group, "filters", 1, &
dims, t % filters, hdf5_err)
end if
! Write universe_bins if present
if (t % n_filter_bins(FILTER_UNIVERSE) > 0) then
dims(1) = t % n_filter_bins(FILTER_UNIVERSE)
call h5ltmake_dataset_int_f(temp_group, "universe_bins", 1, &
dims, t % universe_bins(:) % scalar, hdf5_err)
end if
! Write material_bins if present
if (t % n_filter_bins(FILTER_MATERIAL) > 0) then
dims(1) = t % n_filter_bins(FILTER_MATERIAL)
call h5ltmake_dataset_int_f(temp_group, "material_bins", 1, &
dims, t % material_bins(:) % scalar, hdf5_err)
end if
! Write cell_bins if present
if (t % n_filter_bins(FILTER_CELL) > 0) then
dims(1) = t % n_filter_bins(FILTER_CELL)
call h5ltmake_dataset_int_f(temp_group, "cell_bins", 1, &
dims, t % cell_bins(:) % scalar, hdf5_err)
end if
! Write cellborn_bins if present
if (t % n_filter_bins(FILTER_CELLBORN) > 0) then
dims(1) = t % n_filter_bins(FILTER_CELLBORN)
call h5ltmake_dataset_int_f(temp_group, "cellborn_bins", 1, &
dims, t % cellborn_bins(:) % scalar, hdf5_err)
end if
! Write surface_bins if present
if (t % n_filter_bins(FILTER_SURFACE) > 0) then
dims(1) = t % n_filter_bins(FILTER_SURFACE)
call h5ltmake_dataset_int_f(temp_group, "surface_bins", 1, &
dims, t % surface_bins(:) % scalar, hdf5_err)
end if
! Write incoming energy filter
if (t % n_filter_bins(FILTER_ENERGYIN) > 0) then
dims(1) = t % n_filter_bins(FILTER_ENERGYIN)
call h5ltmake_dataset_double_f(temp_group, "energy_in", 1, &
dims, t % energy_in, hdf5_err)
end if
! Write outgoing energy filter
if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then
dims(1) = t % n_filter_bins(FILTER_ENERGYOUT)
call h5ltmake_dataset_double_f(temp_group, "energy_out", 1, &
dims, t % energy_out, hdf5_err)
end if
! Write mesh information
if (t % n_filter_bins(FILTER_MESH) > 0) then
m => meshes(t % mesh)
dims(1) = m % n_dimension
! Write mesh dimensions
call h5ltmake_dataset_int_f(temp_group, "mesh_dimensions", 1, &
dims, m % dimension, hdf5_err)
! Write mesh lower-left corner
call h5ltmake_dataset_double_f(temp_group, "mesh_lower_left", 1, &
dims, m % origin, hdf5_err)
! Write mesh element width
call h5ltmake_dataset_double_f(temp_group, "mesh_element_width", 1, &
dims, m % width, hdf5_err)
end if
! =======================================================================
! WRITE INFORMATION ON TALLY SCORES
! Create string type of length 12
call h5tcopy_f(H5T_C_S1, string_type, hdf5_err)
call h5tset_size_f(string_type, size_string, hdf5_err)
! Create dataspace and dataset for writing nuclides
dims(1) = t % n_score_bins
call h5screate_simple_f(1, dims, dspace, hdf5_err)
call h5dcreate_f(temp_group, "scores", string_type, dspace, &
dset, hdf5_err)
! Create subspace for writing one element
dims(1) = 1
call h5screate_simple_f(1, dims, subspace, hdf5_err)
! Write list of nuclides
num_elements = 1
do j = 1, t % n_score_bins
! Determine string for this scoring bin
select case (t % score_bins(j) % scalar)
case (SCORE_FLUX)
string = 'flux'
case (SCORE_TOTAL)
string = 'total'
case (SCORE_SCATTER)
string = 'scatter'
case (SCORE_NU_SCATTER)
string = 'nu-scatter'
case (SCORE_SCATTER_1)
string = '1st moment'
case (SCORE_SCATTER_2)
string = '2nd moment'
case (SCORE_SCATTER_3)
string = '3rd moment'
case (SCORE_N_1N)
string = '(n,1n)'
case (SCORE_N_2N)
string = '(n,2n)'
case (SCORE_N_3N)
string = '(n,3n)'
case (SCORE_N_4N)
string = '(n,4n)'
case (SCORE_ABSORPTION)
string = 'absorption'
case (SCORE_FISSION)
string = 'fission'
case (SCORE_NU_FISSION)
string = 'nu-fission'
case (SCORE_CURRENT)
string = 'current'
end select
! Select block within array
coord(1,1) = j
call h5sselect_elements_f(dspace, H5S_SELECT_SET_F, 1, &
num_elements, coord, hdf5_err)
! Write data to partial array
call h5dwrite_f(dset, string_type, string, dims, hdf5_err, &
subspace, dspace)
end do
! Close dataspace and dataset for nuclides
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(subspace, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
! Close group for i-th material
call h5gclose_f(temp_group, hdf5_err)
end do
! Close tallies group
call h5gclose_f(tallies_group, hdf5_err)
end subroutine hdf5_write_tallies
!===============================================================================
! HDF5_WRITE_TIMING
!===============================================================================

22
src/initialize.F90
View file

@ -459,8 +459,8 @@ contains
! =======================================================================
! ADJUST CELL INDICES FOR EACH TALLY
if (t % n_bins(T_CELL) > 0) then
do j = 1, t % n_bins(T_CELL)
if (t % n_filter_bins(FILTER_CELL) > 0) then
do j = 1, t % n_filter_bins(FILTER_CELL)
id = t % cell_bins(j) % scalar
if (dict_has_key(cell_dict, id)) then
t % cell_bins(j) % scalar = dict_get_key(cell_dict, id)
@ -475,8 +475,8 @@ contains
! =======================================================================
! ADJUST SURFACE INDICES FOR EACH TALLY
if (t % n_bins(T_SURFACE) > 0) then
do j = 1, t % n_bins(T_SURFACE)
if (t % n_filter_bins(FILTER_SURFACE) > 0) then
do j = 1, t % n_filter_bins(FILTER_SURFACE)
id = t % surface_bins(j) % scalar
if (dict_has_key(surface_dict, id)) then
t % surface_bins(j) % scalar = dict_get_key(surface_dict, id)
@ -491,8 +491,8 @@ contains
! =======================================================================
! ADJUST UNIVERSE INDICES FOR EACH TALLY
if (t % n_bins(T_UNIVERSE) > 0) then
do j = 1, t % n_bins(T_UNIVERSE)
if (t % n_filter_bins(FILTER_UNIVERSE) > 0) then
do j = 1, t % n_filter_bins(FILTER_UNIVERSE)
id = t % universe_bins(j) % scalar
if (dict_has_key(universe_dict, id)) then
t % universe_bins(j) % scalar = dict_get_key(universe_dict, id)
@ -507,8 +507,8 @@ contains
! =======================================================================
! ADJUST MATERIAL INDICES FOR EACH TALLY
if (t % n_bins(T_MATERIAL) > 0) then
do j = 1, t % n_bins(T_MATERIAL)
if (t % n_filter_bins(FILTER_MATERIAL) > 0) then
do j = 1, t % n_filter_bins(FILTER_MATERIAL)
id = t % material_bins(j) % scalar
if (dict_has_key(material_dict, id)) then
t % material_bins(j) % scalar = dict_get_key(material_dict, id)
@ -523,8 +523,8 @@ contains
! =======================================================================
! ADJUST CELLBORN INDICES FOR EACH TALLY
if (t % n_bins(T_CELLBORN) > 0) then
do j = 1, t % n_bins(T_CELLBORN)
if (t % n_filter_bins(FILTER_CELLBORN) > 0) then
do j = 1, t % n_filter_bins(FILTER_CELLBORN)
id = t % cellborn_bins(j) % scalar
if (dict_has_key(cell_dict, id)) then
t % cellborn_bins(j) % scalar = dict_get_key(cell_dict, id)
@ -539,7 +539,7 @@ contains
! =======================================================================
! ADJUST MESH INDICES FOR EACH TALLY
if (t % n_bins(T_MESH) > 0) then
if (t % n_filter_bins(FILTER_MESH) > 0) then
id = t % mesh
if (dict_has_key(mesh_dict, id)) then
t % mesh = dict_get_key(mesh_dict, id)

165
src/input_xml.F90
View file

@ -617,6 +617,8 @@ contains
integer :: i_mesh ! index in meshes array
integer :: n ! size of arrays in mesh specification
integer :: n_words ! number of words read
integer :: n_filters ! number of filters
integer :: filters(N_FILTER_TYPES) ! temporary list of filters
logical :: file_exists ! does tallies.xml file exist?
character(MAX_LINE_LEN) :: filename
character(MAX_WORD_LEN) :: word
@ -730,13 +732,28 @@ contains
t => tallies(i)
! Allocate arrays for number of bins and stride in scores array
allocate(t % n_bins(TALLY_TYPES))
allocate(t % stride(TALLY_TYPES))
allocate(t % n_filter_bins(N_FILTER_TYPES))
allocate(t % stride(N_FILTER_TYPES))
! Initialize number of bins and stride
t % n_bins = 0
t % n_filter_bins = 0
t % stride = 0
! Initialize filters
n_filters = 0
filters = 0
! Set tally type to volume by default
t % type = TALLY_VOLUME
! It's desirable to use a track-length esimator for tallies since
! generally more events will score to the tally, reducing the
! variance. However, for tallies that require information on
! post-collision parameters (e.g. tally with an energyout filter) the
! analog esimator must be used.
t % estimator = ESTIMATOR_TRACKLENGTH
! Copy material id
t % id = tally_(i) % id
@ -757,7 +774,10 @@ contains
do j = 1, n_words
t % cell_bins(j) % scalar = int(str_to_int(words(j)),4)
end do
t % n_bins(T_CELL) = n_words
t % n_filter_bins(FILTER_CELL) = n_words
n_filters = n_filters + 1
filters(n_filters) = FILTER_CELL
end if
! Read surface filter bins
@ -767,7 +787,10 @@ contains
do j = 1, n_words
t % surface_bins(j) % scalar = int(str_to_int(words(j)),4)
end do
t % n_bins(T_SURFACE) = n_words
t % n_filter_bins(FILTER_SURFACE) = n_words
n_filters = n_filters + 1
filters(n_filters) = FILTER_SURFACE
end if
! Read universe filter bins
@ -777,7 +800,10 @@ contains
do j = 1, n_words
t % universe_bins(j) % scalar = int(str_to_int(words(j)),4)
end do
t % n_bins(T_UNIVERSE) = n_words
t % n_filter_bins(FILTER_UNIVERSE) = n_words
n_filters = n_filters + 1
filters(n_filters) = FILTER_UNIVERSE
end if
! Read material filter bins
@ -787,7 +813,10 @@ contains
do j = 1, n_words
t % material_bins(j) % scalar = int(str_to_int(words(j)),4)
end do
t % n_bins(T_MATERIAL) = n_words
t % n_filter_bins(FILTER_MATERIAL) = n_words
n_filters = n_filters + 1
filters(n_filters) = FILTER_MATERIAL
end if
! Read mesh filter bins
@ -804,7 +833,10 @@ contains
call fatal_error()
end if
t % n_bins(T_MESH) = t % n_bins(T_MESH) + product(m % dimension)
t % n_filter_bins(FILTER_MESH) = t % n_filter_bins(FILTER_MESH) + product(m % dimension)
n_filters = n_filters + 1
filters(n_filters) = FILTER_MESH
end if
! Read birth region filter bins
@ -814,7 +846,10 @@ contains
do j = 1, n_words
t % cellborn_bins(j) % scalar = int(str_to_int(words(j)),4)
end do
t % n_bins(T_CELLBORN) = n_words
t % n_filter_bins(FILTER_CELLBORN) = n_words
n_filters = n_filters + 1
filters(n_filters) = FILTER_CELLBORN
end if
! Read incoming energy filter bins
@ -824,7 +859,10 @@ contains
do j = 1, n_words
t % energy_in(j) = str_to_real(words(j))
end do
t % n_bins(T_ENERGYIN) = n_words - 1
t % n_filter_bins(FILTER_ENERGYIN) = n_words - 1
n_filters = n_filters + 1
filters(n_filters) = FILTER_ENERGYIN
end if
! Read outgoing energy filter bins
@ -834,88 +872,117 @@ contains
do j = 1, n_words
t % energy_out(j) = str_to_real(words(j))
end do
t % n_bins(T_ENERGYOUT) = n_words - 1
t % n_filter_bins(FILTER_ENERGYOUT) = n_words - 1
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
n_filters = n_filters + 1
filters(n_filters) = FILTER_ENERGYOUT
end if
! Allocate and set filters
t % n_filters = n_filters
allocate(t % filters(n_filters))
t % filters = filters(1:n_filters)
! Read macro reactions
if (len_trim(tally_(i) % macros) > 0) then
call split_string(tally_(i) % macros, words, n_words)
allocate(t % macro_bins(n_words))
allocate(t % score_bins(n_words))
do j = 1, n_words
word = words(j)
call lower_case(word)
select case (trim(word))
case ('flux')
t % macro_bins(j) % scalar = MACRO_FLUX
if (t % n_bins(T_ENERGYOUT) > 0) then
t % score_bins(j) % scalar = SCORE_FLUX
if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then
message = "Cannot tally flux with an outgoing energy filter."
call fatal_error()
end if
case ('total')
t % macro_bins(j) % scalar = MACRO_TOTAL
if (t % n_bins(T_ENERGYOUT) > 0) then
t % score_bins(j) % scalar = SCORE_TOTAL
if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then
message = "Cannot tally total reaction rate with an &
&outgoing energy filter."
call fatal_error()
end if
case ('scatter')
t % macro_bins(j) % scalar = MACRO_SCATTER
t % score_bins(j) % scalar = SCORE_SCATTER
case ('nu-scatter')
t % macro_bins(j) % scalar = MACRO_NU_SCATTER
t % score_bins(j) % scalar = SCORE_NU_SCATTER
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
case ('scatter-1')
t % macro_bins(j) % scalar = MACRO_SCATTER_1
t % score_bins(j) % scalar = SCORE_SCATTER_1
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
case ('scatter-2')
t % macro_bins(j) % scalar = MACRO_SCATTER_2
t % score_bins(j) % scalar = SCORE_SCATTER_2
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
case ('scatter-3')
t % macro_bins(j) % scalar = MACRO_SCATTER_3
t % score_bins(j) % scalar = SCORE_SCATTER_3
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
case ('n1n')
t % macro_bins(j) % scalar = MACRO_N_1N
t % score_bins(j) % scalar = SCORE_N_1N
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
case ('n2n')
t % macro_bins(j) % scalar = MACRO_N_2N
t % score_bins(j) % scalar = SCORE_N_2N
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
case ('n3n')
t % macro_bins(j) % scalar = MACRO_N_3N
t % score_bins(j) % scalar = SCORE_N_3N
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
case ('n4n')
t % macro_bins(j) % scalar = MACRO_N_4N
t % score_bins(j) % scalar = SCORE_N_4N
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
case ('absorption')
t % macro_bins(j) % scalar = MACRO_ABSORPTION
if (t % n_bins(T_ENERGYOUT) > 0) then
t % score_bins(j) % scalar = SCORE_ABSORPTION
if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then
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
t % score_bins(j) % scalar = SCORE_FISSION
if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then
message = "Cannot tally fission rate with an outgoing &
&energy filter."
call fatal_error()
end if
case ('nu-fission')
t % macro_bins(j) % scalar = MACRO_NU_FISSION
t % score_bins(j) % scalar = SCORE_NU_FISSION
case ('current')
t % macro_bins(j) % scalar = MACRO_CURRENT
t % surface_current = .true.
t % score_bins(j) % scalar = SCORE_CURRENT
t % type = TALLY_SURFACE_CURRENT
! 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 scoring functions in the same &
&tally as surface currents. Separate other scoring &
&functions into a distinct tally."
call fatal_error()
end if
! Check to make sure that only the mesh filter was specified
!!$ if (t % mesh == 0 .or. t % n_bins(T_MESH) /= &
!!$ product(t % n_bins, t % n_bins > 0)) then
!!$ message = "Surface currents must be used with a mesh filter only."
!!$ call fatal_error()
!!$ end if
! Since the number of bins for the mesh filter was already set
! assuming it was a flux tally, we need to adjust the number of
! bins
t % n_bins(T_MESH) = t % n_bins(T_MESH) - product(m % dimension)
t % n_filter_bins(FILTER_MESH) = t % n_filter_bins(FILTER_MESH) &
- product(m % dimension)
! Get pointer to mesh
id = t % mesh
@ -925,19 +992,19 @@ contains
! We need to increase the dimension by one since we also need
! currents coming into and out of the boundary mesh cells.
if (size(m % dimension) == 2) then
t % n_bins(T_MESH) = t % n_bins(T_MESH) + &
product(m % dimension + 1) * 4
t % n_filter_bins(FILTER_MESH) = t % n_filter_bins(FILTER_MESH) &
+ product(m % dimension + 1) * 4
elseif (size(m % dimension) == 3) then
t % n_bins(T_MESH) = t % n_bins(T_MESH) + &
product(m % dimension + 1) * 6
t % n_filter_bins(FILTER_MESH) = t % n_filter_bins(FILTER_MESH) &
+ product(m % dimension + 1) * 6
end if
case default
message = "Unknown macro reaction: " // trim(words(j))
message = "Unknown scoring function: " // trim(words(j))
call fatal_error()
end select
end do
t % n_macro_bins = n_words
t % n_score_bins = n_words
end if
end do

2
src/intercycle.F90
View file

@ -299,7 +299,7 @@ contains
t => tallies(i)
n = t % n_total_bins
m = t % n_macro_bins
m = t % n_score_bins
n_bins = n*m
allocate(tally_temp(n,m))

48
src/output.F90
View file

@ -575,9 +575,9 @@ contains
write(unit_,*) 'Tally ' // to_str(t % id)
if (t % n_bins(T_CELL) > 0) then
if (t % n_filter_bins(FILTER_CELL) > 0) then
string = ""
do i = 1, t % n_bins(T_CELL)
do i = 1, t % n_filter_bins(FILTER_CELL)
id = t % cell_bins(i) % scalar
c => cells(id)
string = trim(string) // ' ' // trim(to_str(c % id))
@ -585,9 +585,9 @@ contains
write(unit_, *) ' Cell Bins:' // trim(string)
end if
if (t % n_bins(T_SURFACE) > 0) then
if (t % n_filter_bins(FILTER_SURFACE) > 0) then
string = ""
do i = 1, t % n_bins(T_SURFACE)
do i = 1, t % n_filter_bins(FILTER_SURFACE)
id = t % surface_bins(i) % scalar
s => surfaces(id)
string = trim(string) // ' ' // trim(to_str(s % id))
@ -595,9 +595,9 @@ contains
write(unit_, *) ' Surface Bins:' // trim(string)
end if
if (t % n_bins(T_UNIVERSE) > 0) then
if (t % n_filter_bins(FILTER_UNIVERSE) > 0) then
string = ""
do i = 1, t % n_bins(T_UNIVERSE)
do i = 1, t % n_filter_bins(FILTER_UNIVERSE)
id = t % universe_bins(i) % scalar
u => universes(id)
string = trim(string) // ' ' // trim(to_str(u % id))
@ -605,9 +605,9 @@ contains
write(unit_, *) ' Material Bins:' // trim(string)
end if
if (t % n_bins(T_MATERIAL) > 0) then
if (t % n_filter_bins(FILTER_MATERIAL) > 0) then
string = ""
do i = 1, t % n_bins(T_MATERIAL)
do i = 1, t % n_filter_bins(FILTER_MATERIAL)
id = t % material_bins(i) % scalar
m => materials(id)
string = trim(string) // ' ' // trim(to_str(m % id))
@ -615,7 +615,7 @@ contains
write(unit_, *) ' Material Bins:' // trim(string)
end if
if (t % n_bins(T_MESH) > 0) then
if (t % n_filter_bins(FILTER_MESH) > 0) then
string = ""
id = t % mesh
sm => meshes(id)
@ -626,9 +626,9 @@ contains
write(unit_, *) ' Mesh Bins:' // trim(string)
end if
if (t % n_bins(T_CELLBORN) > 0) then
if (t % n_filter_bins(FILTER_CELLBORN) > 0) then
string = ""
do i = 1, t % n_bins(T_CELLBORN)
do i = 1, t % n_filter_bins(FILTER_CELLBORN)
id = t % cellborn_bins(i) % scalar
c => cells(id)
string = trim(string) // ' ' // trim(to_str(c % id))
@ -636,39 +636,39 @@ contains
write(unit_, *) ' Birth Region Bins:' // trim(string)
end if
if (t % n_bins(T_ENERGYIN) > 0) then
if (t % n_filter_bins(FILTER_ENERGYIN) > 0) then
string = ""
do i = 1, t % n_bins(T_ENERGYIN) + 1
do i = 1, t % n_filter_bins(FILTER_ENERGYIN) + 1
string = trim(string) // ' ' // trim(to_str(&
t % energy_in(i)))
end do
write(unit_,*) ' Incoming Energy Bins:' // trim(string)
end if
if (t % n_bins(T_ENERGYOUT) > 0) then
if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then
string = ""
do i = 1, t % n_bins(T_ENERGYOUT) + 1
do i = 1, t % n_filter_bins(FILTER_ENERGYOUT) + 1
string = trim(string) // ' ' // trim(to_str(&
t % energy_out(i)))
end do
write(unit_,*) ' Outgoing Energy Bins:' // trim(string)
end if
if (t % n_macro_bins > 0) then
if (t % n_score_bins > 0) then
string = ""
do i = 1, t % n_macro_bins
select case (t % macro_bins(i) % scalar)
case (MACRO_FLUX)
do i = 1, t % n_score_bins
select case (t % score_bins(i) % scalar)
case (SCORE_FLUX)
string = trim(string) // ' flux'
case (MACRO_TOTAL)
case (SCORE_TOTAL)
string = trim(string) // ' total'
case (MACRO_SCATTER)
case (SCORE_SCATTER)
string = trim(string) // ' scatter'
case (MACRO_ABSORPTION)
case (SCORE_ABSORPTION)
string = trim(string) // ' absorption'
case (MACRO_FISSION)
case (SCORE_FISSION)
string = trim(string) // ' fission'
case (MACRO_NU_FISSION)
case (SCORE_NU_FISSION)
string = trim(string) // ' nu-fission'
end select
end do

8
src/particle_header.F90
View file

@ -1,6 +1,6 @@
module particle_header
use constants, only: NEUTRON, ONE, NONE
use constants, only: NEUTRON, ONE, NONE, ZERO
use geometry_header, only: BASE_UNIVERSE
implicit none
@ -23,9 +23,8 @@ module particle_header
real(8) :: xyz(3)
real(8) :: uvw(3)
! Pointers to next (lower) and previous (higher) universe
! Pointer to next (more local) set of coordinates
type(LocalCoord), pointer :: next => null()
type(LocalCoord), pointer :: prev => null()
end type LocalCoord
!===============================================================================
@ -53,6 +52,9 @@ module particle_header
real(8) :: last_wgt ! last particle weight
real(8) :: last_E ! last energy
! What event last took place
integer :: event
! Post-collision physical data
integer :: n_bank ! number of fission sites banked

60
src/physics.F90
View file

@ -16,7 +16,8 @@ module physics
use random_lcg, only: prn
use search, only: binary_search
use string, only: to_str
use tally, only: score_tally, score_surface_current
use tally, only: score_analog_tally, score_tracklength_tally, &
score_surface_current
implicit none
@ -45,7 +46,6 @@ contains
! Particle couldn't be located
if (.not. found_cell) then
message = "Could not locate particle " // trim(to_str(p % id))
print *, p % coord0 % xyz
call fatal_error()
end if
@ -67,7 +67,6 @@ contains
! Initialize number of events to zero
n_event = 0
! find energy index, interpolation factor
do while (p % alive)
! Calculate microscopic and macroscopic cross sections -- note: if the
@ -96,15 +95,22 @@ contains
coord => coord % next
end do
! Score track-length tallies
if (tallies_on) call score_tracklength_tally(p, distance)
if (d_collision > d_boundary) then
last_cell = p % coord % cell
p % coord % cell = NONE
if (lattice_crossed /= NONE) then
! Particle crosses lattice boundary
p % surface = NONE
call cross_lattice(p, lattice_crossed)
p % event = EVENT_LATTICE
else
! Particle crosses surface
p % surface = surface_crossed
call cross_surface(p, last_cell)
p % event = EVENT_SURFACE
end if
else
! collision
@ -149,10 +155,6 @@ contains
type(Particle), pointer :: p
integer :: MT ! ENDF reaction number
logical :: scattered ! was this a scattering reaction?
logical :: fissioned ! was this a fission reaction?
! Store pre-collision particle properties
p % last_wgt = p % wgt
p % last_E = p % E
@ -167,13 +169,13 @@ contains
if (tallies_on) call score_surface_current(p)
! Sample nuclide/reaction for the material the particle is in
call sample_reaction(p, MT)
call sample_reaction(p)
! Display information about collision
if (verbosity >= 10 .or. trace) then
message = " " // trim(reaction_name(MT)) // ". Energy = " // &
trim(to_str(p % E * 1e6_8)) // " eV."
call write_message()
!!$ message = " " // trim(reaction_name(MT)) // ". Energy = " // &
!!$ trim(to_str(p % E * 1e6_8)) // " eV."
!!$ call write_message()
end if
! check for very low energy
@ -183,21 +185,12 @@ contains
call warning()
end if
! Check if particle scattered or fissioned
if (survival_biasing) then
fissioned = .false.
scattered = .true.
else
fissioned = is_fission(MT)
scattered = is_scatter(MT)
end if
! Score collision estimator tallies for any macro tallies -- this is done
! after a collision has occurred rather than before because we need
! information on the outgoing energy for any tallies with an outgoing energy
! filter
if (tallies_on) call score_tally(p, scattered, fissioned)
if (tallies_on) call score_analog_tally(p)
! Reset number of particles banked during collision
p % n_bank = 0
@ -212,10 +205,9 @@ contains
! disappearance are treated implicitly.
!===============================================================================
subroutine sample_reaction(p, MT)
subroutine sample_reaction(p)
type(Particle), pointer :: p
integer, intent(out) :: MT ! ENDF MT number of reaction that occured
integer :: i ! index over nuclides in a material
integer :: index_nuclide ! index in nuclides array
@ -282,7 +274,7 @@ contains
! See if disappearance reaction happens
if (prob > cutoff) then
p % alive = .false.
MT = N_DISAPPEAR
p % event = EVENT_ABSORB
return
end if
end if
@ -338,16 +330,24 @@ contains
if (micro_xs(index_nuclide) % use_ptable) then
! In the case that the particle is in the unresolved resonance
! region and probability tables are being used, we should only
! check the first fission reaction
prob = prob + micro_xs(index_nuclide) % fission
if (prob > cutoff) then
rxn => nuc % reactions(nuc % index_fission(1))
call create_fission_sites(p, index_nuclide, rxn, .true.)
p % alive = .false.
MT = rxn % MT
p % event = EVENT_FISSION
return
end if
else
! With no probability tables, we need to loop through each fission
! reaction type
do i = 1, nuc % n_fission
rxn => nuc % reactions(nuc % index_fission(i))
@ -366,7 +366,7 @@ contains
if (prob > cutoff) then
call create_fission_sites(p, index_nuclide, rxn, .true.)
p % alive = .false.
MT = rxn % MT
p % event = EVENT_FISSION
return
end if
end do
@ -419,8 +419,6 @@ contains
end if
! Set MT to be returned
MT = 2
else
! =======================================================================
! INELASTIC SCATTERING
@ -461,10 +459,12 @@ contains
! Perform collision physics for inelastics scattering
call inelastic_scatter(p, nuc, rxn)
! Set MT to be returned
MT = rxn % MT
end if
! If we made it this far, it means that a scattering reaction took place
! since the absorption and fission blocks had return statements.
p % event = EVENT_SCATTER
end subroutine sample_reaction
!===============================================================================

932
src/tally.F90
View file

File diff suppressed because it is too large Load diff

35
src/tally_header.F90
View file

@ -58,12 +58,17 @@ module tally_header
type TallyObject
! Basic data
integer :: id
integer :: type
real(8) :: volume
logical :: surface_current = .false.
integer :: id ! user-defined identifier
integer :: type ! volume, surface current
integer :: estimator ! collision, track-length
real(8) :: volume ! volume of region
! Tally bin specifications
! Information about what filters should be used
integer :: n_filters
integer, allocatable :: filters(:)
! Filter bin specifications
type(TallyFilter), pointer :: universe_bins(:) => null()
type(TallyFilter), pointer :: material_bins(:) => null()
@ -74,24 +79,24 @@ module tally_header
real(8), allocatable :: energy_in(:)
real(8), allocatable :: energy_out(:)
! Number of bins for each filter
integer :: n_total_bins = 0
! Total number of filter bins
integer :: n_total_bins = 0
! The following attributes do not necessarily need to be stored but they
! greatly simplify logic in many places. n_bins gives the number of bins
! for each filter type, e.g. n_bins(T_CELL) would be the size of
! cell_bins. The stride attribute is used for determining the index in the
! scores array for a bin combination. Since multiple dimensions are mapped
! onto one dimension in the scores array, the stride attribute gives the
! stride for a given filter type within the scores array
! for each filter type, e.g. n_filter_bins(FILTER_CELL) would be the size
! of cell_bins. The stride attribute is used for determining the index in
! the scores array for a bin combination. Since multiple dimensions are
! mapped onto one dimension in the scores array, the stride attribute gives
! the stride for a given filter type within the scores array
integer, allocatable :: n_bins(:)
integer, allocatable :: n_filter_bins(:)
integer, allocatable :: stride(:)
! Macroscopic properties to score
type(TallyFilter), pointer :: macro_bins(:) => null()
integer :: n_macro_bins = 0
type(TallyFilter), pointer :: score_bins(:) => null()
integer :: n_score_bins = 0
! Scores for each bin -- the most natural way to have scores would be to
! have a dimension for each different type of bin, but older Fortran

2
src/utils/analyze_source.py
View file

@ -34,6 +34,8 @@ for sourceFile in source:
commentChar = '!'
elif ending == 'f90':
commentChar = '!'
elif ending == 'F90':
commentChar = '!'
for line in open(sourceFile, 'r'):
line = line.strip()