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

Merging with upstream.

Browse source
This commit is contained in:
Adam Nelson 2012-12-05 10:13:03 -05:00
commit eb7cc7d895
29 changed files with 771 additions and 734 deletions
Download patch file Download diff file Expand all files Collapse all files

5
docs/source/releasenotes/notes_0.5.0.rst
View file

@ -32,6 +32,9 @@ New Features
Bug Fixes
---------
- 737b90_: Coincident surfaces from separate universes / particle traveling
tangent to a surface.
- a819b4_: Output of surface neighbors in summary.out file.
- b11696_: Reading long attribute lists in XML input.
- 2bd46a_: Search for tallying nuclides when no default_xs specified.
- 7a1f08_: Fix word wrapping when writing messages.
@ -39,6 +42,8 @@ Bug Fixes
- 6f8d9d_: Set default tally labels.
- 6a3a5e_: Fix problem with corner-crossing in lattices.
.. _737b90: https://github.com/mit-crpg/openmc/commit/737b90
.. _a819b4: https://github.com/mit-crpg/openmc/commit/a819b4
.. _b11696: https://github.com/mit-crpg/openmc/commit/b11696
.. _2bd46a: https://github.com/mit-crpg/openmc/commit/2bd46a
.. _7a1f08: https://github.com/mit-crpg/openmc/commit/7a1f08

62
docs/source/usersguide/input.rst
View file

@ -69,37 +69,6 @@ deviation.
*Default*: off
``<criticality>`` Element
-------------------------
The ``<criticality>`` element indicates that a criticality calculation should be
performed. It has the following attributes/sub-elements:
:batches:
The total number of batches, where each batch corresponds to multiple
fission source iterations. Batching is done to eliminate correlation between
realizations of random variables.
*Default*: None
:generations_per_batch:
The number of total fission source iterations per batch.
*Default*: 1
:inactive:
The number of inactive batches. In general, the starting cycles in a
criticality calculation can not be used to contribute to tallies since the
fission source distribution and eigenvalue are generally not converged
immediately.
*Default*: None
:particles:
The number of neutrons to simulate per fission source iteration.
*Default*: None
.. _cross_sections:
``<cross_sections>`` Element
@ -130,6 +99,37 @@ default. This element has the following attributes/sub-elements:
*Default*: 1.0
``<eigenvalue>`` Element
------------------------
The ``<eigenvalue>`` element indicates that a :math:`k`-eigenvalue calculation
should be performed. It has the following attributes/sub-elements:
:batches:
The total number of batches, where each batch corresponds to multiple
fission source iterations. Batching is done to eliminate correlation between
realizations of random variables.
*Default*: None
:generations_per_batch:
The number of total fission source iterations per batch.
*Default*: 1
:inactive:
The number of inactive batches. In general, the starting cycles in a
criticality calculation can not be used to contribute to tallies since the
fission source distribution and eigenvalue are generally not converged
immediately.
*Default*: None
:particles:
The number of neutrons to simulate per fission source iteration.
*Default*: None
``<energy_grid>`` Element
-------------------------

36
src/DEPENDENCIES
View file

@ -90,23 +90,23 @@ cmfd_snes_solver.o: constants.o
cmfd_snes_solver.o: global.o
cmfd_snes_solver.o: string.o
criticality.o: cmfd_execute.o
criticality.o: constants.o
criticality.o: error.o
criticality.o: global.o
criticality.o: hdf5_interface.o
criticality.o: math.o
criticality.o: mesh.o
criticality.o: mesh_header.o
criticality.o: output.o
criticality.o: physics.o
criticality.o: random_lcg.o
criticality.o: search.o
criticality.o: source.o
criticality.o: state_point.o
criticality.o: string.o
criticality.o: tally.o
criticality.o: timing.o
eigenvalue.o: cmfd_execute.o
eigenvalue.o: constants.o
eigenvalue.o: error.o
eigenvalue.o: global.o
eigenvalue.o: hdf5_interface.o
eigenvalue.o: math.o
eigenvalue.o: mesh.o
eigenvalue.o: mesh_header.o
eigenvalue.o: output.o
eigenvalue.o: physics.o
eigenvalue.o: random_lcg.o
eigenvalue.o: search.o
eigenvalue.o: source.o
eigenvalue.o: state_point.o
eigenvalue.o: string.o
eigenvalue.o: tally.o
eigenvalue.o: timing.o
cross_section.o: ace_header.o
cross_section.o: constants.o
@ -243,7 +243,7 @@ interpolation.o: search.o
interpolation.o: string.o
main.o: constants.o
main.o: criticality.o
main.o: eigenvalue.o
main.o: finalize.o
main.o: fixed_source.o
main.o: global.o

2
src/OBJECTS
View file

@ -14,10 +14,10 @@ cmfd_power_solver.o \
cmfd_prod_operator.o \
cmfd_snes_solver.o \
cross_section.o \
criticality.o \
datatypes.o \
datatypes_header.o \
doppler.o \
eigenvalue.o \
endf.o \
endf_header.o \
energy_grid.o \

2
src/bank_header.F90
View file

@ -3,7 +3,7 @@ module bank_header
implicit none
!===============================================================================
! BANK is used for storing fission sites in criticality calculations. Since all
! BANK is used for storing fission sites in eigenvalue calculations. Since all
! the state information of a neutron is not needed, this type allows sites to be
! stored with less memory
!===============================================================================

40
src/cmfd_data.F90
View file

@ -61,7 +61,7 @@ contains
meshes
use mesh, only: mesh_indices_to_bin
use mesh_header, only: StructuredMesh
use tally_header, only: TallyObject, TallyScore
use tally_header, only: TallyObject
integer :: nx ! number of mesh cells in x direction
integer :: ny ! number of mesh cells in y direction
@ -155,7 +155,7 @@ contains
score_index = sum((t % matching_bins - 1) * t%stride) + 1
! get flux
flux = t % scores(1,score_index) % sum
flux = t % results(1,score_index) % sum
cmfd % flux(h,i,j,k) = flux
! detect zero flux
@ -166,13 +166,13 @@ contains
end if
! get total rr and convert to total xs
cmfd % totalxs(h,i,j,k) = t % scores(2,score_index) % sum / flux
cmfd % totalxs(h,i,j,k) = t % results(2,score_index) % sum / flux
! get p1 scatter rr and convert to p1 scatter xs
cmfd % p1scattxs(h,i,j,k) = t % scores(3,score_index) % sum / flux
cmfd % p1scattxs(h,i,j,k) = t % results(3,score_index) % sum / flux
! extract diffusion coefficient tally
cmfd % diffusion(h,i,j,k) = t % scores(4,score_index) % sum / flux
cmfd % diffusion(h,i,j,k) = t % results(4,score_index) % sum / flux
! calculate diffusion coefficient
! cmfd % diffcof(h,i,j,k) = ONE/(3.0_8*cmfd%totalxs(h,i,j,k))
@ -206,16 +206,16 @@ contains
score_index = sum((t % matching_bins - 1) * t%stride) + 1
! get scattering
cmfd % scattxs(h,g,i,j,k) = t % scores(1,score_index) % sum /&
cmfd % scattxs(h,g,i,j,k) = t % results(1,score_index) % sum /&
cmfd % flux(h,i,j,k)
! get nu-fission
cmfd % nfissxs(h,g,i,j,k) = t % scores(2,score_index) % sum /&
cmfd % nfissxs(h,g,i,j,k) = t % results(2,score_index) % sum /&
cmfd % flux(h,i,j,k)
! bank source
cmfd % openmc_src(g,i,j,k) = cmfd % openmc_src(g,i,j,k) + &
t % scores(2,score_index) % sum
t % results(2,score_index) % sum
end do INGROUP
@ -232,60 +232,60 @@ contains
(/ i-1, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_RIGHT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
cmfd % current(1,h,i,j,k) = t % scores(1,score_index) % sum
cmfd % current(1,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_RIGHT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
cmfd % current(2,h,i,j,k) = t % scores(1,score_index) % sum
cmfd % current(2,h,i,j,k) = t % results(1,score_index) % sum
! right surface
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_RIGHT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
cmfd % current(3,h,i,j,k) = t % scores(1,score_index) % sum
cmfd % current(3,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_RIGHT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
cmfd % current(4,h,i,j,k) = t % scores(1,score_index) % sum
cmfd % current(4,h,i,j,k) = t % results(1,score_index) % sum
! back surface
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j-1, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_FRONT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
cmfd % current(5,h,i,j,k) = t % scores(1,score_index) % sum
cmfd % current(5,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_FRONT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
cmfd % current(6,h,i,j,k) = t % scores(1,score_index) % sum
cmfd % current(6,h,i,j,k) = t % results(1,score_index) % sum
! front surface
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_FRONT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
cmfd % current(7,h,i,j,k) = t % scores(1,score_index) % sum
cmfd % current(7,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_FRONT
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
cmfd % current(8,h,i,j,k) = t % scores(1,score_index) % sum
cmfd % current(8,h,i,j,k) = t % results(1,score_index) % sum
! bottom surface
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k-1 /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_TOP
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
cmfd % current(9,h,i,j,k) = t % scores(1,score_index) % sum
cmfd % current(9,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_TOP
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
cmfd % current(10,h,i,j,k) = t % scores(1,score_index) % sum
cmfd % current(10,h,i,j,k) = t % results(1,score_index) % sum
! top surface
t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_TOP
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming
cmfd % current(11,h,i,j,k) = t % scores(1,score_index) % sum
cmfd % current(11,h,i,j,k) = t % results(1,score_index) % sum
t % matching_bins(i_filter_surf) = OUT_TOP
score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing
cmfd % current(12,h,i,j,k) = t % scores(1,score_index) % sum
cmfd % current(12,h,i,j,k) = t % results(1,score_index) % sum
end if TALLY

4
src/cmfd_execute.F90
View file

@ -478,7 +478,7 @@ contains
use global, only: n_user_tallies, n_tallies, tallies, message
use output, only: write_message
use tally, only: reset_score
use tally, only: reset_result
integer :: i ! loop counter
@ -491,7 +491,7 @@ contains
! reset that tally
tallies(i) % n_realizations = 0
call reset_score(tallies(i) % scores)
call reset_result(tallies(i) % results)
end do

2
src/cmfd_input.F90
View file

@ -173,7 +173,7 @@ contains
use mesh_header, only: StructuredMesh
use string
use tally, only: setup_active_cmfdtallies
use tally_header, only: TallyObject, TallyScore, TallyFilter
use tally_header, only: TallyObject, TallyFilter
use xml_data_cmfd_t
integer :: i ! loop counter

3
src/constants.F90
View file

@ -32,6 +32,7 @@ module constants
! User for precision in geometry
real(8), parameter :: FP_PRECISION = 1e-14_8
real(8), parameter :: FP_REL_PRECISION = 1e-5_8
real(8), parameter :: FP_COINCIDENT = 1e-12_8
! Maximum number of collisions/crossings
integer, parameter :: MAX_EVENTS = 10000
@ -348,7 +349,7 @@ module constants
! Running modes
integer, parameter :: &
MODE_FIXEDSOURCE = 1, & ! Fixed source mode
MODE_CRITICALITY = 2, & ! Criticality mode
MODE_EIGENVALUE = 2, & ! K eigenvalue mode
MODE_PLOTTING = 3, & ! Plotting mode
MODE_TALLIES = 4 ! Tally results mode

52
src/criticality.F90 → src/eigenvalue.F90
View file

@ -1,4 +1,4 @@
module criticality
module eigenvalue
#ifdef MPI
use mpi
@ -26,18 +26,21 @@ module criticality
use hdf5_interface, only: hdf5_write_state_point
#endif
private
public :: run_eigenvalue
contains
!===============================================================================
! RUN_CRITICALITY encompasses all the main logic where iterations are performed
! over the batches, generations, and histories.
! RUN_EIGENVALUE encompasses all the main logic where iterations are performed
! over the batches, generations, and histories in a k-eigenvalue calculation.
!===============================================================================
subroutine run_criticality()
subroutine run_eigenvalue()
integer(8) :: i ! index over individual particles
if (master) call header("CRITICALITY TRANSPORT SIMULATION", level=1)
if (master) call header("K EIGENVALUE SIMULATION", level=1)
! Allocate particle
allocate(p)
@ -45,6 +48,9 @@ contains
! Display column titles
call print_columns()
! Turn on inactive timer
call timer_start(time_inactive)
! ==========================================================================
! LOOP OVER BATCHES
BATCH_LOOP: do current_batch = 1, n_batches
@ -99,7 +105,7 @@ contains
if (master) call header("SIMULATION FINISHED", level=1)
end subroutine run_criticality
end subroutine run_eigenvalue
!===============================================================================
! INITIALIZE_BATCH
@ -113,19 +119,22 @@ contains
! Reset total starting particle weight used for normalizing tallies
total_weight = ZERO
if (current_batch == n_inactive + 1) then
! Switch from inactive batch timer to active batch timer
call timer_stop(time_inactive)
call timer_start(time_active)
! Enable active batches (and tallies_on if it hasn't been enabled)
active_batches = .true.
tallies_on = .true.
! Add user tallies to active tallies list
call setup_active_usertallies()
end if
! check CMFD initialize batch
if (cmfd_run) call cmfd_init_batch()
if (current_batch == n_inactive + 1) then
! This will start the active timer at the first non-inactive batch
! (including batch 1 if there are no inactive batches).
call timer_start(time_active)
elseif (current_batch == 1) then
! If there are inactive batches, start the inactive timer on the first
! batch.
call timer_start(time_inactive)
end if
end subroutine initialize_batch
!===============================================================================
@ -184,15 +193,6 @@ contains
end if
end do
! Turn tallies on once inactive batches are complete
if (current_batch == n_inactive) then
tallies_on = .true.
active_batches = .true.
call timer_stop(time_inactive)
call timer_start(time_active)
call setup_active_usertallies()
end if
end subroutine finalize_batch
!===============================================================================
@ -692,4 +692,4 @@ contains
end subroutine count_source_for_ufs
end module criticality
end module eigenvalue

57
src/geometry.F90
View file

@ -8,7 +8,7 @@ module geometry
use output, only: write_message
use particle_header, only: LocalCoord, deallocate_coord
use string, only: to_str
use tally, only: score_surface_current, add_to_score
use tally, only: score_surface_current
implicit none
@ -47,7 +47,7 @@ contains
! Determine the specified sense of the surface in the cell and the actual
! sense of the particle with respect to the surface
s => surfaces(abs(i_surface))
actual_sense = sense(s, p % coord % xyz)
actual_sense = sense(s)
specified_sense = (c % surfaces(i) > 0)
! Compare sense of point to specified sense
@ -272,10 +272,9 @@ contains
call score_surface_current()
end if
if (tallies_on) then
! Score to global leakage tally
call add_to_score(global_tallies(LEAKAGE), p % wgt)
end if
! Score to global leakage tally
if (tallies_on) global_tallies(LEAKAGE) % value = &
global_tallies(LEAKAGE) % value + p % wgt
! Display message
if (verbosity >= 10 .or. trace) then
@ -481,14 +480,32 @@ contains
! ==========================================================================
! COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS
! Remove lower coordinate levels and assignment of surface
p % surface = NONE
p % coord => p % coord0
call deallocate_coord(p % coord % next)
call find_cell(found)
! Couldn't find next cell anywhere!
if ((.not. found) .and. (run_mode /= MODE_PLOTTING)) then
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()
if (run_mode /= MODE_PLOTTING .and. (.not. found)) then
! If a cell is still not found, there are two possible causes: 1) there is
! a void in the model, and 2) the particle hit a surface at a tangent. If
! the particle is really traveling tangent to a surface, if we move it
! forward a tiny bit it should fix the problem.
p % coord => p % coord0
call deallocate_coord(p % coord % next)
p % coord % xyz = p % coord % xyz + TINY_BIT * p % coord % uvw
call find_cell(found)
! Couldn't find next cell anywhere! This probably means there is an actual
! undefined region in the geometry.
if (.not. found) then
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
end if
end subroutine cross_surface
@ -1142,10 +1159,9 @@ contains
! is in.
!===============================================================================
function sense(surf, xyz) result(s)
recursive function sense(surf) result(s)
type(Surface), pointer :: surf ! surface
real(8), intent(in) :: xyz(3) ! coordinates of particle
logical :: s ! sense of particle
real(8) :: x,y,z ! coordinates of particle
@ -1164,9 +1180,9 @@ contains
real(8) :: r ! radius for quadratic surfaces
real(8) :: x1,y1,z1 ! upper-right corner of box
x = xyz(1)
y = xyz(2)
z = xyz(3)
x = p % coord % xyz(1)
y = p % coord % xyz(2)
z = p % coord % xyz(3)
select case (surf % type)
case (SURF_PX)
@ -1320,7 +1336,12 @@ contains
end select
! Check which side of surface the point is on
if (func > 0) then
if (abs(func) < FP_COINCIDENT) then
! Particle may be coincident with this surface. Artifically move the
! particle forward a tiny bit.
p % coord % xyz = p % coord % xyz + TINY_BIT * p % coord % uvw
s = sense(surf)
elseif (func > 0) then
s = .true.
else
s = .false.

20
src/global.F90
View file

@ -13,7 +13,7 @@ module global
use particle_header, only: Particle
use plot_header, only: Plot
use source_header, only: ExtSource
use tally_header, only: TallyObject, TallyMap, TallyScore
use tally_header, only: TallyObject, TallyMap, TallyResult
use timing, only: Timer
#ifdef MPI
@ -111,7 +111,7 @@ module global
! 3) track-length estimate of k-eff
! 4) leakage fraction
type(TallyScore), target :: global_tallies(N_GLOBAL_TALLIES)
type(TallyResult), target :: global_tallies(N_GLOBAL_TALLIES)
! Tally map structure
type(TallyMap), allocatable :: tally_maps(:)
@ -154,7 +154,7 @@ module global
type(ListInt), pointer :: active_tallies => null()
! ============================================================================
! CRITICALITY SIMULATION VARIABLES
! EIGENVALUE SIMULATION VARIABLES
integer(8) :: n_particles = 0 ! # of particles per generation
integer :: n_batches ! # of batches
@ -208,7 +208,7 @@ module global
logical :: mpi_enabled = .false. ! is MPI in use and initialized?
integer :: mpi_err ! MPI error code
integer :: MPI_BANK ! MPI datatype for fission bank
integer :: MPI_TALLYSCORE ! MPI datatype for TallyScore
integer :: MPI_TALLYRESULT ! MPI datatype for TallyResult
! No reduction at end of batch
logical :: reduce_tallies = .true.
@ -240,17 +240,17 @@ module global
! HDF5 VARIABLES
#ifdef HDF5
integer(HID_T) :: hdf5_output_file ! identifier for output file
integer(HID_T) :: hdf5_tallyscore_t ! Compound type for TallyScore
integer(HID_T) :: hdf5_bank_t ! Compound type for Bank
integer(HID_T) :: hdf5_integer8_t ! type for integer(8)
integer :: hdf5_err ! error flag
integer(HID_T) :: hdf5_output_file ! identifier for output file
integer(HID_T) :: hdf5_tallyresult_t ! Compound type for TallyResult
integer(HID_T) :: hdf5_bank_t ! Compound type for Bank
integer(HID_T) :: hdf5_integer8_t ! type for integer(8)
integer :: hdf5_err ! error flag
#endif
! ============================================================================
! MISCELLANEOUS VARIABLES
! Mode to run in (fixed source, criticality, plotting, etc)
! Mode to run in (fixed source, eigenvalue, plotting, etc)
integer :: run_mode = NONE
! Restart run

60
src/hdf5_interface.F90
View file

@ -30,20 +30,20 @@ contains
subroutine hdf5_initialize()
type(TallyScore), target :: tmp(2) ! temporary TallyScore
type(Bank), target :: tmpb(2) ! temporary Bank
integer(HID_T) :: coordinates_t ! HDF5 type for 3 reals
integer(HSIZE_T) :: dims(1) = (/3/) ! size of coordinates
type(TallyResult), target :: tmp(2) ! temporary TallyResult
type(Bank), target :: tmpb(2) ! temporary Bank
integer(HID_T) :: coordinates_t ! HDF5 type for 3 reals
integer(HSIZE_T) :: dims(1) = (/3/) ! size of coordinates
! Initialize FORTRAN interface.
call h5open_f(hdf5_err)
! Create the compound datatype for TallyScore
! Create the compound datatype for TallyResult
call h5tcreate_f(H5T_COMPOUND_F, h5offsetof(c_loc(tmp(1)), &
c_loc(tmp(2))), hdf5_tallyscore_t, hdf5_err)
call h5tinsert_f(hdf5_tallyscore_t, "sum", h5offsetof(c_loc(tmp(1)), &
c_loc(tmp(2))), hdf5_tallyresult_t, hdf5_err)
call h5tinsert_f(hdf5_tallyresult_t, "sum", h5offsetof(c_loc(tmp(1)), &
c_loc(tmp(1)%sum)), H5T_NATIVE_DOUBLE, hdf5_err)
call h5tinsert_f(hdf5_tallyscore_t, "sum_sq", h5offsetof(c_loc(tmp(1)), &
call h5tinsert_f(hdf5_tallyresult_t, "sum_sq", h5offsetof(c_loc(tmp(1)), &
c_loc(tmp(1)%sum_sq)), H5T_NATIVE_DOUBLE, hdf5_err)
! Create compound type for xyz and uvw
@ -73,7 +73,7 @@ contains
subroutine hdf5_finalize()
! Release compound datatypes
call h5tclose_f(hdf5_tallyscore_t, hdf5_err)
call h5tclose_f(hdf5_tallyresult_t, hdf5_err)
! Close FORTRAN interface.
call h5close_f(hdf5_err)
@ -94,8 +94,8 @@ contains
! Write header information
call hdf5_write_header()
! Write criticality information
if (run_mode == MODE_CRITICALITY) then
! Write eigenvalue information
if (run_mode == MODE_EIGENVALUE) then
! Need to write integer(8)'s using double instead since there is no H5LT
! call for making a dataset of type long
call hdf5_write_double(hdf5_output_file, "n_particles", real(n_particles,8))
@ -865,8 +865,8 @@ contains
! Write out current batch number
call hdf5_write_integer(hdf5_state_point, "current_batch", current_batch)
! Write out information for criticality run
if (run_mode == MODE_CRITICALITY) then
! Write out information for eigenvalue run
if (run_mode == MODE_EIGENVALUE) then
call hdf5_write_integer(hdf5_state_point, "n_inactive", n_inactive)
call hdf5_write_integer(hdf5_state_point, "gen_per_batch", gen_per_batch)
@ -1031,10 +1031,10 @@ contains
! Write global tallies
dims(1) = N_GLOBAL_TALLIES
call h5screate_simple_f(1, dims, dspace, hdf5_err)
call h5dcreate_f(hdf5_state_point, "global_tallies", hdf5_tallyscore_t, &
call h5dcreate_f(hdf5_state_point, "global_tallies", hdf5_tallyresult_t, &
dspace, dset, hdf5_err)
f_ptr = c_loc(global_tallies(1))
CALL h5dwrite_f(dset, hdf5_tallyscore_t, f_ptr, hdf5_err)
CALL h5dwrite_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err)
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
@ -1043,7 +1043,7 @@ contains
call hdf5_write_integer(tallies_group, "tallies_present", 1)
! Write tally sum and sum_sq
TALLY_SCORES: do i = 1, n_tallies
TALLY_RESULTS: do i = 1, n_tallies
! Get pointer to tally
t => tallies(i)
@ -1052,18 +1052,18 @@ contains
temp_group, hdf5_err)
! Write sum and sum_sq for each bin
dims2 = shape(t % scores)
dims2 = shape(t % results)
call h5screate_simple_f(2, dims2, dspace, hdf5_err)
call h5dcreate_f(temp_group, "values", hdf5_tallyscore_t, &
call h5dcreate_f(temp_group, "results", hdf5_tallyresult_t, &
dspace, dset, hdf5_err)
f_ptr = c_loc(t % scores(1, 1))
CALL h5dwrite_f(dset, hdf5_tallyscore_t, f_ptr, hdf5_err)
f_ptr = c_loc(t % results(1, 1))
CALL h5dwrite_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err)
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
! Close group for the i-th tally
call h5gclose_f(temp_group, hdf5_err)
end do TALLY_SCORES
end do TALLY_RESULTS
else
! Indicate that tallies are off
call hdf5_write_integer(tallies_group, "tallies_present", 0)
@ -1132,8 +1132,8 @@ contains
! Read batch number to restart at
call hdf5_read_integer(hdf5_state_point, "current_batch", restart_batch)
! Read information specific to criticality run
if (mode == MODE_CRITICALITY) then
! Read information specific to eigenvalue run
if (mode == MODE_EIGENVALUE) then
call hdf5_read_integer(hdf5_state_point, "n_inactive", n_inactive)
call hdf5_read_integer(hdf5_state_point, "gen_per_batch", gen_per_batch)
@ -1153,7 +1153,7 @@ contains
! Read global tallies
f_ptr = c_loc(global_tallies(1))
call h5dread_f(dset, hdf5_tallyscore_t, f_ptr, hdf5_err)
call h5dread_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err)
! Close global tallies dataset
call h5dclose_f(dset, hdf5_err)
@ -1167,14 +1167,14 @@ contains
call hdf5_read_integer(tally_group, "n_realizations", &
tallies(i) % n_realizations)
! Open dataset for tally values
call h5dopen_f(tally_group, "values", dset, hdf5_err)
! Open dataset for tally results
call h5dopen_f(tally_group, "results", dset, hdf5_err)
! Read sum and sum_sq for each tally bin
f_ptr = c_loc(tallies(i) % scores(1,1))
call h5dread_f(dset, hdf5_tallyscore_t, f_ptr, hdf5_err)
f_ptr = c_loc(tallies(i) % results(1,1))
call h5dread_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err)
! Close dataset for tally values
! Close dataset for tally results
call h5dclose_f(dset, hdf5_err)
! Close tally group
@ -1191,7 +1191,7 @@ contains
f_ptr = c_loc(source_bank(1))
call h5dread_f(dset, hdf5_bank_t, f_ptr, hdf5_err)
! Close dataset for tally values
! Close dataset for source bank
call h5dclose_f(dset, hdf5_err)
! Close HDF5 state point file

38
src/initialize.F90
View file

@ -120,7 +120,7 @@ contains
! Allocate banks and create source particles -- for a fixed source
! calculation, the external source distribution is sampled during the
! run, not at initialization
if (run_mode == MODE_CRITICALITY) then
if (run_mode == MODE_EIGENVALUE) then
call allocate_banks()
if (.not. restart_run) call initialize_source()
end if
@ -169,14 +169,14 @@ contains
integer :: bank_types(4) ! Datatypes
integer(MPI_ADDRESS_KIND) :: bank_disp(4) ! Displacements
integer :: temp_type ! temporary derived type
integer :: score_blocks(1) ! Count for each datatype
integer :: score_types(1) ! Datatypes
integer(MPI_ADDRESS_KIND) :: score_disp(1) ! Displacements
integer(MPI_ADDRESS_KIND) :: score_base_disp ! Base displacement
integer(MPI_ADDRESS_KIND) :: lower_bound ! Lower bound for TallyScore
integer(MPI_ADDRESS_KIND) :: extent ! Extent for TallyScore
integer :: result_blocks(1) ! Count for each datatype
integer :: result_types(1) ! Datatypes
integer(MPI_ADDRESS_KIND) :: result_disp(1) ! Displacements
integer(MPI_ADDRESS_KIND) :: result_base_disp ! Base displacement
integer(MPI_ADDRESS_KIND) :: lower_bound ! Lower bound for TallyResult
integer(MPI_ADDRESS_KIND) :: extent ! Extent for TallyResult
type(Bank) :: b
type(TallyScore) :: ts
type(TallyResult) :: tr
! Indicate that MPI is turned on
mpi_enabled = .true.
@ -215,29 +215,29 @@ contains
call MPI_TYPE_COMMIT(MPI_BANK, mpi_err)
! ==========================================================================
! CREATE MPI_TALLYSCORE TYPE
! CREATE MPI_TALLYRESULT TYPE
! Determine displacements for MPI_BANK type
call MPI_GET_ADDRESS(ts % value, score_base_disp, mpi_err)
call MPI_GET_ADDRESS(ts % sum, score_disp(1), mpi_err)
call MPI_GET_ADDRESS(tr % value, result_base_disp, mpi_err)
call MPI_GET_ADDRESS(tr % sum, result_disp(1), mpi_err)
! Adjust displacements
score_disp = score_disp - score_base_disp
result_disp = result_disp - result_base_disp
! Define temporary type for tallyscore
score_blocks = (/ 2 /)
score_types = (/ MPI_REAL8 /)
call MPI_TYPE_CREATE_STRUCT(1, score_blocks, score_disp, score_types, &
! Define temporary type for TallyResult
result_blocks = (/ 2 /)
result_types = (/ MPI_REAL8 /)
call MPI_TYPE_CREATE_STRUCT(1, result_blocks, result_disp, result_types, &
temp_type, mpi_err)
! Adjust lower-bound and extent of type for tally score
lower_bound = 0
extent = score_disp(1) + 16
extent = result_disp(1) + 16
call MPI_TYPE_CREATE_RESIZED(temp_type, lower_bound, extent, &
MPI_TALLYSCORE, mpi_err)
MPI_TALLYRESULT, mpi_err)
! Commit derived type for tally scores
call MPI_TYPE_COMMIT(MPI_TALLYSCORE, mpi_err)
call MPI_TYPE_COMMIT(MPI_TALLYRESULT, mpi_err)
end subroutine initialize_mpi
#endif

335
src/input_xml.F90
View file

@ -19,7 +19,7 @@ module input_xml
implicit none
type(DictionaryII), pointer :: & ! used to count how many cells each
cells_in_univ_dict => null() ! universe contains
cells_in_univ_dict => null() ! universe contains
contains
@ -92,7 +92,7 @@ contains
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."
&or in the CROSS_SECTIONS environment variable."
call fatal_error()
else
path_cross_sections = trim(env_variable)
@ -102,32 +102,41 @@ contains
end if
! Make sure that either criticality or fixed source was specified
if (criticality_ % batches == 0 .and. fixed_source_ % batches == 0) then
message = "Number of batches on <criticality> or <fixed_source> " &
// "tag was zero."
if (eigenvalue_ % batches == 0 .and. fixed_source_ % batches == 0 &
.and. criticality_ % batches == 0) then
message = "Number of batches on <eigenvalue> or <fixed_source> &
&tag was zero."
call fatal_error()
end if
! Criticality information
! Check for old <criticality> tag
if (criticality_ % batches > 0) then
eigenvalue_ = criticality_
message = "The <criticality> element has been deprecated and &
&replaced by <eigenvalue>."
call warning()
end if
! Eigenvalue information
if (eigenvalue_ % batches > 0) then
! Set run mode
if (run_mode == NONE) run_mode = MODE_CRITICALITY
if (run_mode == NONE) run_mode = MODE_EIGENVALUE
! Check number of particles
if (len_trim(criticality_ % particles) == 0) then
if (len_trim(eigenvalue_ % particles) == 0) then
message = "Need to specify number of particles per cycles."
call fatal_error()
end if
! If the number of particles was specified as a command-line argument, we
! don't set it here
if (n_particles == 0) n_particles = str_to_int(criticality_ % particles)
if (n_particles == 0) n_particles = str_to_int(eigenvalue_ % particles)
! Copy batch and generation information
n_batches = criticality_ % batches
n_inactive = criticality_ % inactive
n_batches = eigenvalue_ % batches
n_inactive = eigenvalue_ % inactive
n_active = n_batches - n_inactive
gen_per_batch = criticality_ % generations_per_batch
gen_per_batch = eigenvalue_ % generations_per_batch
! Allocate array for batch keff and entropy
allocate(k_batch(n_batches))
@ -168,9 +177,6 @@ contains
call fatal_error()
end if
! Turn on tallies if no inactive batches
if (n_inactive == 0) tallies_on = .true.
! Copy random number seed if specified
if (seed_ > 0) seed = seed_
@ -202,7 +208,7 @@ contains
inquire(FILE=path_source, EXIST=file_exists)
if (.not. file_exists) then
message = "Binary source file '" // trim(path_source) // &
"' does not exist!"
"' does not exist!"
call fatal_error()
end if
@ -221,7 +227,7 @@ contains
coeffs_reqd = 3
case default
message = "Invalid spatial distribution for external source: " &
// trim(source_ % space % type)
// trim(source_ % space % type)
call fatal_error()
end select
@ -235,11 +241,11 @@ contains
! Read parameters for spatial distribution
if (n < coeffs_reqd) then
message = "Not enough parameters specified for spatial " &
// "distribution of external source."
// "distribution of external source."
call fatal_error()
elseif (n > coeffs_reqd) then
message = "Too many parameters specified for spatial " &
// "distribution of external source."
// "distribution of external source."
call fatal_error()
elseif (n > 0) then
allocate(external_source % params_space(n))
@ -266,7 +272,7 @@ contains
external_source % type_angle = SRC_ANGLE_TABULAR
case default
message = "Invalid angular distribution for external source: " &
// trim(source_ % angle % type)
// trim(source_ % angle % type)
call fatal_error()
end select
@ -280,11 +286,11 @@ contains
! Read parameters for angle distribution
if (n < coeffs_reqd) then
message = "Not enough parameters specified for angle " &
// "distribution of external source."
// "distribution of external source."
call fatal_error()
elseif (n > coeffs_reqd) then
message = "Too many parameters specified for angle " &
// "distribution of external source."
// "distribution of external source."
call fatal_error()
elseif (n > 0) then
allocate(external_source % params_angle(n))
@ -314,7 +320,7 @@ contains
external_source % type_energy = SRC_ENERGY_TABULAR
case default
message = "Invalid energy distribution for external source: " &
// trim(source_ % energy % type)
// trim(source_ % energy % type)
call fatal_error()
end select
@ -328,11 +334,11 @@ contains
! Read parameters for energy distribution
if (n < coeffs_reqd) then
message = "Not enough parameters specified for energy " &
// "distribution of external source."
// "distribution of external source."
call fatal_error()
elseif (n > coeffs_reqd) then
message = "Too many parameters specified for energy " &
// "distribution of external source."
// "distribution of external source."
call fatal_error()
elseif (n > 0) then
allocate(external_source % params_energy(n))
@ -370,10 +376,10 @@ contains
! Check to make sure enough values were supplied
if (size(entropy_(1) % lower_left) /= 3) then
message = "Need to specify (x,y,z) coordinates of lower-left corner &
&of Shannon entropy mesh."
&of Shannon entropy mesh."
elseif (size(entropy_(1) % upper_right) /= 3) then
message = "Need to specify (x,y,z) coordinates of upper-right corner &
&of Shannon entropy mesh."
&of Shannon entropy mesh."
end if
! Allocate mesh object and coordinates on mesh
@ -388,7 +394,7 @@ contains
! Check on values provided
if (.not. all(entropy_mesh % upper_right > entropy_mesh % lower_left)) then
message = "Upper-right coordinate must be greater than lower-left &
&coordinate for Shannon entropy mesh."
&coordinate for Shannon entropy mesh."
call fatal_error()
end if
@ -399,7 +405,7 @@ contains
! If so, make sure proper number of values were given
if (size(entropy_(1) % dimension) /= 3) then
message = "Dimension of entropy mesh must be given as three &
&integers."
&integers."
call fatal_error()
end if
@ -420,13 +426,13 @@ contains
! Check to make sure enough values were supplied
if (size(uniform_fs_(1) % lower_left) /= 3) then
message = "Need to specify (x,y,z) coordinates of lower-left corner &
&of UFS mesh."
&of UFS mesh."
elseif (size(uniform_fs_(1) % upper_right) /= 3) then
message = "Need to specify (x,y,z) coordinates of upper-right corner &
&of UFS mesh."
&of UFS mesh."
elseif (size(uniform_fs_(1) % dimension) /= 3) then
message = "Dimension of UFS mesh must be given as three &
&integers."
&integers."
call fatal_error()
end if
@ -450,13 +456,13 @@ contains
! Check on values provided
if (.not. all(ufs_mesh % upper_right > ufs_mesh % lower_left)) then
message = "Upper-right coordinate must be greater than lower-left &
&coordinate for UFS mesh."
&coordinate for UFS mesh."
call fatal_error()
end if
! Calculate width
ufs_mesh % width = (ufs_mesh % upper_right - &
ufs_mesh % lower_left) / ufs_mesh % dimension
ufs_mesh % lower_left) / ufs_mesh % dimension
! Calculate volume fraction of each cell
ufs_mesh % volume_frac = ONE/real(product(ufs_mesh % dimension),8)
@ -466,7 +472,7 @@ contains
! Allocate source_frac
allocate(source_frac(1, ufs_mesh % dimension(1), &
ufs_mesh % dimension(2), ufs_mesh % dimension(3)))
ufs_mesh % dimension(2), ufs_mesh % dimension(3)))
end if
! Check if the user has specified to write state points
@ -488,7 +494,7 @@ contains
n_state_points = n_batches / state_point_(1) % interval
allocate(statepoint_batch(n_state_points))
statepoint_batch = (/ (state_point_(1) % interval * i, i = 1, &
n_state_points) /)
n_state_points) /)
else
! If neither were specified, write state point at last batch
n_state_points = 1
@ -498,7 +504,7 @@ contains
! Check if the user has specified to write binary source file
if (trim(state_point_(1) % source_separate) == 'on') &
source_separate = .true.
source_separate = .true.
else
! If no <state_point> tag was present, by default write state point at
! last batch only
@ -629,7 +635,7 @@ contains
! Check to make sure that either material or fill was specified
if (c % material == NONE .and. c % fill == NONE) then
message = "Neither material nor fill was specified for cell " // &
trim(to_str(c % id))
trim(to_str(c % id))
call fatal_error()
end if
@ -643,7 +649,7 @@ contains
! Check to make sure that surfaces were specified
if (.not. associated(cell_(i) % surfaces)) then
message = "No surfaces specified for cell " // &
trim(to_str(c % id))
trim(to_str(c % id))
call fatal_error()
end if
@ -659,7 +665,7 @@ contains
! another universe
if (c % fill == NONE) then
message = "Cannot apply a rotation to cell " // trim(to_str(&
c % id)) // " because it is not filled with another universe"
c % id)) // " because it is not filled with another universe"
call fatal_error()
end if
@ -667,7 +673,7 @@ contains
n = size(cell_(i) % rotation)
if (n /= 3) then
message = "Incorrect number of rotation parameters on cell " // &
to_str(c % id)
to_str(c % id)
call fatal_error()
end if
@ -679,13 +685,13 @@ contains
! Calculate rotation matrix based on angles given
allocate(c % rotation(3,3))
c % rotation = reshape((/ &
cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta), &
-cos(phi)*sin(psi) + sin(phi)*sin(theta)*cos(psi), &
cos(phi)*cos(psi) + sin(phi)*sin(theta)*sin(psi), &
sin(phi)*cos(theta), &
sin(phi)*sin(psi) + cos(phi)*sin(theta)*cos(psi), &
-sin(phi)*cos(psi) + cos(phi)*sin(theta)*sin(psi), &
cos(phi)*cos(theta) /), (/ 3,3 /))
cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta), &
-cos(phi)*sin(psi) + sin(phi)*sin(theta)*cos(psi), &
cos(phi)*cos(psi) + sin(phi)*sin(theta)*sin(psi), &
sin(phi)*cos(theta), &
sin(phi)*sin(psi) + cos(phi)*sin(theta)*cos(psi), &
-sin(phi)*cos(psi) + cos(phi)*sin(theta)*sin(psi), &
cos(phi)*cos(theta) /), (/ 3,3 /))
end if
! Translation vector
@ -694,7 +700,7 @@ contains
! another universe
if (c % fill == NONE) then
message = "Cannot apply a translation to cell " // trim(to_str(&
c % id)) // " because it is not filled with another universe"
c % id)) // " because it is not filled with another universe"
call fatal_error()
end if
@ -702,7 +708,7 @@ contains
n = size(cell_(i) % translation)
if (n /= 3) then
message = "Incorrect number of translation parameters on cell " &
// to_str(c % id)
// to_str(c % id)
call fatal_error()
end if
@ -814,11 +820,11 @@ contains
n = size(surface_(i) % coeffs)
if (n < coeffs_reqd) then
message = "Not enough coefficients specified for surface: " // &
trim(to_str(s % id))
trim(to_str(s % id))
call fatal_error()
elseif (n > coeffs_reqd) then
message = "Too many coefficients specified for surface: " // &
trim(to_str(s % id))
trim(to_str(s % id))
call fatal_error()
else
allocate(s % coeffs(n))
@ -839,7 +845,7 @@ contains
s % bc = BC_PERIODIC
case default
message = "Unknown boundary condition '" // trim(word) // &
"' specified on surface " // trim(to_str(s % id))
"' specified on surface " // trim(to_str(s % id))
call fatal_error()
end select
@ -888,7 +894,7 @@ contains
! Read lattice lower-left location
if (size(lattice_(i) % dimension) /= size(lattice_(i) % lower_left)) then
message = "Number of entries on <lower_left> must be the same as " // &
"the number of entries on <dimension>."
"the number of entries on <dimension>."
call fatal_error()
end if
l % x0 = lattice_(i) % lower_left(1)
@ -897,7 +903,7 @@ contains
! Read lattice widths
if (size(lattice_(i) % width) /= size(lattice_(i) % lower_left)) then
message = "Number of entries on <width> must be the same as " // &
"the number of entries on <lower_left>."
"the number of entries on <lower_left>."
call fatal_error()
end if
l % width_x = lattice_(i) % width(1)
@ -998,7 +1004,7 @@ contains
sum_density = .false.
if (val <= ZERO) then
message = "Need to specify a positive density on material " // &
trim(to_str(mat % id)) // "."
trim(to_str(mat % id)) // "."
call fatal_error()
end if
@ -1015,7 +1021,7 @@ contains
mat % density = 1.0e-24 * val
case default
message = "Unkwown units '" // trim(material_(i) % density % units) &
// "' specified on material " // trim(to_str(mat % id))
// "' specified on material " // trim(to_str(mat % id))
call fatal_error()
end select
end if
@ -1026,7 +1032,7 @@ contains
! Check to ensure material has at least one nuclide
if (.not. associated(material_(i) % nuclides)) then
message = "No nuclides specified on material " // &
trim(to_str(mat % id))
trim(to_str(mat % id))
call fatal_error()
end if
@ -1044,7 +1050,7 @@ contains
! Check for empty name on nuclide
if (len_trim(nuc % name) == 0) then
message = "No name specified on nuclide in material " // &
trim(to_str(mat % id))
trim(to_str(mat % id))
call fatal_error()
end if
@ -1052,7 +1058,7 @@ contains
if (len_trim(nuc % xs) == 0) then
if (default_xs == '') then
message = "No cross section specified for nuclide in material " &
// trim(to_str(mat % id))
// trim(to_str(mat % id))
call fatal_error()
else
nuc % xs = default_xs
@ -1066,7 +1072,7 @@ contains
! Check that this nuclide is listed in the cross_sections.xml file
if (.not. dict_has_key(xs_listing_dict, name)) then
message = "Could not find nuclide " // trim(name) // &
" in cross_sections.xml file!"
" in cross_sections.xml file!"
call fatal_error()
end if
@ -1074,7 +1080,7 @@ contains
n = len_trim(name)
if (name(n:n) /= 'c') then
message = "Cross-section table " // trim(name) // &
" is not a continuous-energy neutron table."
" is not a continuous-energy neutron table."
call fatal_error()
end if
@ -1099,11 +1105,11 @@ contains
! weight percents were specified
if (nuc % ao == ZERO .and. nuc % wo == ZERO) then
message = "No atom or weight percent specified for nuclide " // &
trim(name)
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)
&nuclide: " // trim(name)
call fatal_error()
end if
@ -1118,9 +1124,9 @@ contains
! Check to make sure either all atom percents or all weight percents are
! given
if (.not. (all(mat % atom_density > ZERO) .or. &
all(mat % atom_density < ZERO))) then
all(mat % atom_density < ZERO))) then
message = "Cannot mix atom and weight percents in material " // &
to_str(mat % id)
to_str(mat % id)
call fatal_error()
end if
@ -1141,7 +1147,7 @@ contains
! Check that this nuclide is listed in the cross_sections.xml file
if (.not. dict_has_key(xs_listing_dict, name)) then
message = "Could not find S(a,b) table " // trim(name) // &
" in cross_sections.xml file!"
" in cross_sections.xml file!"
call fatal_error()
end if
mat % has_sab_table = .true.
@ -1304,7 +1310,7 @@ contains
! Check that dimensions are all greater than zero
if (any(mesh_(i) % dimension <= 0)) then
message = "All entries on the <dimension> element for a tally mesh &
&must be positive."
&must be positive."
call fatal_error()
end if
@ -1314,24 +1320,24 @@ contains
! Read mesh lower-left corner location
if (m % n_dimension /= size(mesh_(i) % lower_left)) then
message = "Number of entries on <lower_left> must be the same as &
&the number of entries on <dimension>."
&the number of entries on <dimension>."
call fatal_error()
end if
m % lower_left = mesh_(i) % lower_left
! Make sure either upper-right or width was specified
if (associated(mesh_(i) % upper_right) .and. &
associated(mesh_(i) % width)) then
associated(mesh_(i) % width)) then
message = "Cannot specify both <upper_right> and <width> on a &
&tally mesh."
&tally mesh."
call fatal_error()
end if
! Make sure either upper-right or width was specified
if (.not. associated(mesh_(i) % upper_right) .and. &
.not. associated(mesh_(i) % width)) then
.not. associated(mesh_(i) % width)) then
message = "Must specify either <upper_right> and <width> on a &
&tally mesh."
&tally mesh."
call fatal_error()
end if
@ -1339,7 +1345,7 @@ contains
! Check to ensure width has same dimensions
if (size(mesh_(i) % width) /= size(mesh_(i) % lower_left)) then
message = "Number of entries on <width> must be the same as the &
&number of entries on <lower_left>."
&number of entries on <lower_left>."
call fatal_error()
end if
@ -1357,14 +1363,14 @@ contains
! Check to ensure width has same dimensions
if (size(mesh_(i) % upper_right) /= size(mesh_(i) % lower_left)) then
message = "Number of entries on <upper_right> must be the same as &
&the number of entries on <lower_left>."
&the number of entries on <lower_left>."
call fatal_error()
end if
! Check that upper-right is above lower-left
if (any(mesh_(i) % upper_right < mesh_(i) % lower_left)) then
message = "The <upper_right> coordinates must be greater than the &
&<lower_left> coordinates on a tally mesh."
&<lower_left> coordinates on a tally mesh."
call fatal_error()
end if
@ -1413,7 +1419,7 @@ contains
if (size(tally_(i) % filters) > 0) then
message = "Tally filters should be specified with multiple <filter> &
&elements. Did you forget to change your <filters> element?"
&elements. Did you forget to change your <filters> element?"
call fatal_error()
end if
@ -1445,7 +1451,7 @@ contains
allocate(t % filters(j) % int_bins(n_words))
do k = 1, n_words
t % filters(j) % int_bins(k) = int(str_to_int(&
tally_(i) % filter(j) % bins(k)),4)
tally_(i) % filter(j) % bins(k)),4)
end do
case ('cellborn')
@ -1459,7 +1465,7 @@ contains
allocate(t % filters(j) % int_bins(n_words))
do k = 1, n_words
t % filters(j) % int_bins(k) = int(str_to_int(&
tally_(i) % filter(j) % bins(k)),4)
tally_(i) % filter(j) % bins(k)),4)
end do
case ('material')
@ -1473,7 +1479,7 @@ contains
allocate(t % filters(j) % int_bins(n_words))
do k = 1, n_words
t % filters(j) % int_bins(k) = int(str_to_int(&
tally_(i) % filter(j) % bins(k)),4)
tally_(i) % filter(j) % bins(k)),4)
end do
case ('universe')
@ -1487,7 +1493,7 @@ contains
allocate(t % filters(j) % int_bins(n_words))
do k = 1, n_words
t % filters(j) % int_bins(k) = int(str_to_int(&
tally_(i) % filter(j) % bins(k)),4)
tally_(i) % filter(j) % bins(k)),4)
end do
case ('surface')
@ -1501,7 +1507,7 @@ contains
allocate(t % filters(j) % int_bins(n_words))
do k = 1, n_words
t % filters(j) % int_bins(k) = int(str_to_int(&
tally_(i) % filter(j) % bins(k)),4)
tally_(i) % filter(j) % bins(k)),4)
end do
case ('mesh')
@ -1523,7 +1529,7 @@ contains
m => meshes(i_mesh)
else
message = "Could not find mesh " // trim(to_str(id)) // &
" specified on tally " // trim(to_str(t % id))
" specified on tally " // trim(to_str(t % id))
call fatal_error()
end if
@ -1547,7 +1553,7 @@ contains
allocate(t % filters(j) % real_bins(n_words))
do k = 1, n_words
t % filters(j) % real_bins(k) = str_to_real(&
tally_(i) % filter(j) % bins(k))
tally_(i) % filter(j) % bins(k))
end do
case ('energyout')
@ -1561,7 +1567,7 @@ contains
allocate(t % filters(j) % real_bins(n_words))
do k = 1, n_words
t % filters(j) % real_bins(k) = str_to_real(&
tally_(i) % filter(j) % bins(k))
tally_(i) % filter(j) % bins(k))
end do
end select
@ -1575,9 +1581,9 @@ contains
! Check that both cell and surface weren't specified
if (t % find_filter(FILTER_CELL) > 0 .and. &
t % find_filter(FILTER_SURFACE) > 0) then
t % find_filter(FILTER_SURFACE) > 0) then
message = "Cannot specify both cell and surface filters for tally " &
// trim(to_str(t % id))
// trim(to_str(t % id))
call fatal_error()
end if
@ -1596,7 +1602,7 @@ contains
! Set bins to 1, 2, 3, ..., n_nuclides_total, -1
t % nuclide_bins(1:n_nuclides_total) = &
(/ (j, j=1, n_nuclides_total) /)
(/ (j, j=1, n_nuclides_total) /)
t % nuclide_bins(n_nuclides_total + 1) = -1
! Set number of nuclide bins
@ -1624,7 +1630,7 @@ contains
key_list => dict_keys(nuclide_dict)
do while (associated(key_list))
if (starts_with(key_list % data % key, &
tally_(i) % nuclides(j))) then
tally_(i) % nuclides(j))) then
word = key_list % data % key
exit
end if
@ -1636,8 +1642,8 @@ contains
! Check if no nuclide was found
if (.not. associated(key_list)) then
message = "Could not find the nuclide " // trim(&
tally_(i) % nuclides(j)) // " specified in tally " &
// trim(to_str(t % id)) // " in any material."
tally_(i) % nuclides(j)) // " specified in tally " &
// trim(to_str(t % id)) // " in any material."
call fatal_error()
end if
deallocate(key_list)
@ -1650,7 +1656,7 @@ contains
! Check to make sure nuclide specified is in problem
if (.not. dict_has_key(nuclide_dict, word)) then
message = "The nuclide " // trim(word) // " from tally " // &
trim(to_str(t % id)) // " is not present in any material."
trim(to_str(t % id)) // " is not present in any material."
call fatal_error()
end if
@ -1749,7 +1755,7 @@ contains
case ('flux')
! Prohibit user from tallying flux for an individual nuclide
if (.not. (t % n_nuclide_bins == 1 .and. &
t % nuclide_bins(1) == -1)) then
t % nuclide_bins(1) == -1)) then
message = "Cannot tally flux for an individual nuclide."
call fatal_error()
end if
@ -1763,7 +1769,7 @@ contains
t % score_bins(j) = SCORE_TOTAL
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
message = "Cannot tally total reaction rate with an &
&outgoing energy filter."
&outgoing energy filter."
call fatal_error()
end if
case ('scatter')
@ -1824,14 +1830,14 @@ contains
t % score_bins(j) = SCORE_ABSORPTION
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
message = "Cannot tally absorption rate with an outgoing &
&energy filter."
&energy filter."
call fatal_error()
end if
case ('fission')
t % score_bins(j) = SCORE_FISSION
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
message = "Cannot tally fission rate with an outgoing &
&energy filter."
&energy filter."
call fatal_error()
end if
case ('nu-fission')
@ -1844,8 +1850,8 @@ contains
! for this tally
if (n_words > 1) then
message = "Cannot tally other scoring functions in the same &
&tally as surface currents. Separate other scoring &
&functions into a distinct tally."
&tally as surface currents. Separate other scoring &
&functions into a distinct tally."
call fatal_error()
end if
@ -1877,13 +1883,13 @@ contains
t % filters(t % n_filters) % type = FILTER_SURFACE
t % filters(t % n_filters) % n_bins = 2 * m % n_dimension
allocate(t % filters(t % n_filters) % int_bins(&
2 * m % n_dimension))
2 * m % n_dimension))
if (m % n_dimension == 2) then
t % filters(t % n_filters) % int_bins = (/ IN_RIGHT, &
OUT_RIGHT, IN_FRONT, OUT_FRONT /)
OUT_RIGHT, IN_FRONT, OUT_FRONT /)
elseif (m % n_dimension == 3) then
t % filters(t % n_filters) % int_bins = (/ IN_RIGHT, &
OUT_RIGHT, IN_FRONT, OUT_FRONT, IN_TOP, OUT_TOP /)
OUT_RIGHT, IN_FRONT, OUT_FRONT, IN_TOP, OUT_TOP /)
end if
t % find_filter(FILTER_SURFACE) = t % n_filters
@ -1892,15 +1898,14 @@ contains
case default
message = "Unknown scoring function: " // &
trim(tally_(i) % scores(l))
trim(tally_(i) % scores(j))
call fatal_error()
end select
end do
t % n_score_bins = n_scores
t % n_nonPN_score_bins = n_words
t % n_score_bins = n_words
else
message = "No <scores> specified on tally " // trim(to_str(t % id)) &
// "."
// "."
call fatal_error()
end if
@ -1918,7 +1923,7 @@ contains
! tally needs post-collision information
if (t % estimator == ESTIMATOR_ANALOG) then
message = "Cannot use track-length estimator for tally " &
// to_str(t % id)
// to_str(t % id)
call fatal_error()
end if
@ -1927,7 +1932,7 @@ contains
case default
message = "Invalid estimator '" // trim(tally_(i) % estimator) &
// "' on tally " // to_str(t % id)
// "' on tally " // to_str(t % id)
call fatal_error()
end select
end if
@ -2031,7 +2036,7 @@ contains
! Set output file path
pl % path_plot = trim(path_input) // trim(to_str(pl % id)) // &
"_" // trim(plot_(i) % filename) // ".ppm"
"_" // trim(plot_(i) % filename) // ".ppm"
! Copy plot pixel size
if (size(plot_(i) % pixels) == 2) then
@ -2047,7 +2052,7 @@ contains
pl % not_found % rgb = plot_(i) % background
else
message = "Bad background RGB " &
// "in plot " // trim(to_str(pl % id))
// "in plot " // trim(to_str(pl % id))
call fatal_error()
end if
else
@ -2062,7 +2067,7 @@ contains
! pl % type = PLOT_TYPE_POINTS
case default
message = "Unsupported plot type '" // plot_(i) % type &
// "' in plot " // trim(to_str(pl % id))
// "' in plot " // trim(to_str(pl % id))
call fatal_error()
end select
@ -2076,7 +2081,7 @@ contains
pl % basis = PLOT_BASIS_YZ
case default
message = "Unsupported plot basis '" // plot_(i) % basis &
// "' in plot " // trim(to_str(pl % id))
// "' in plot " // trim(to_str(pl % id))
call fatal_error()
end select
@ -2085,7 +2090,7 @@ contains
pl % origin = plot_(i) % origin
else
message = "Origin must be length 3 " &
// "in plot " // trim(to_str(pl % id))
// "in plot " // trim(to_str(pl % id))
call fatal_error()
end if
@ -2097,7 +2102,7 @@ contains
pl % width(2) = plot_(i) % width(2)
else
message = "Bad plot width " &
// "in plot " // trim(to_str(pl % id))
// "in plot " // trim(to_str(pl % id))
call fatal_error()
end if
@ -2125,7 +2130,7 @@ contains
case default
message = "Unsupported plot color type '" // plot_(i) % color &
// "' in plot " // trim(to_str(pl % id))
// "' in plot " // trim(to_str(pl % id))
call fatal_error()
end select
@ -2135,7 +2140,7 @@ contains
do j = 1, n_cols
if (size(plot_(i) % col_spec_(j) % rgb) /= 3) then
message = "Bad RGB " &
// "in plot " // trim(to_str(pl % id))
// "in plot " // trim(to_str(pl % id))
call fatal_error()
end if
@ -2147,7 +2152,7 @@ contains
pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb
else
message = "Could not find cell " // trim(to_str(col_id)) // &
" specified in plot " // trim(to_str(pl % id))
" specified in plot " // trim(to_str(pl % id))
call fatal_error()
end if
@ -2157,7 +2162,7 @@ contains
pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb
else
message = "Could not find material " // trim(to_str(col_id)) // &
" specified in plot " // trim(to_str(pl % id))
" specified in plot " // trim(to_str(pl % id))
call fatal_error()
end if
@ -2169,7 +2174,7 @@ contains
if (associated(plot_(i) % mask_)) then
if (size(plot_(i) % mask_) > 1) then
message = "Mutliple masks" // &
" specified in plot " // trim(to_str(pl % id))
" specified in plot " // trim(to_str(pl % id))
call fatal_error()
else if (.not. size(plot_(i) % mask_) == 0) then
do j=1,size(pl % colors)
@ -2204,10 +2209,10 @@ contains
! Check if cross_sections.xml exists
inquire(FILE=path_cross_sections, EXIST=file_exists)
if (.not. file_exists) then
! Could not find cross_sections.xml file
message = "Cross sections XML file '" // trim(path_cross_sections) // &
! Could not find cross_sections.xml file
message = "Cross sections XML file '" // trim(path_cross_sections) // &
"' does not exist!"
call fatal_error()
call fatal_error()
end if
message = "Reading cross sections XML file..."
@ -2223,25 +2228,25 @@ contains
call read_xml_file_cross_sections_t(path_cross_sections)
if (len_trim(directory_) > 0) then
! Copy directory information if present
directory = trim(directory_)
! Copy directory information if present
directory = trim(directory_)
else
! If no directory is listed in cross_sections.xml, by default select the
! directory in which the cross_sections.xml file resides
i = index(path_cross_sections, "/", BACK=.true.)
directory = path_cross_sections(1:i)
! If no directory is listed in cross_sections.xml, by default select the
! directory in which the cross_sections.xml file resides
i = index(path_cross_sections, "/", BACK=.true.)
directory = path_cross_sections(1:i)
end if
! determine whether binary/ascii
if (filetype_ == 'ascii') then
filetype = ASCII
filetype = ASCII
elseif (filetype_ == 'binary') then
filetype = BINARY
filetype = BINARY
elseif (len_trim(filetype_) == 0) then
filetype = ASCII
filetype = ASCII
else
message = "Unknown filetype in cross_sections.xml: " // trim(filetype_)
call fatal_error()
message = "Unknown filetype in cross_sections.xml: " // trim(filetype_)
call fatal_error()
end if
! copy default record length and entries for binary files
@ -2250,57 +2255,57 @@ contains
! Allocate xs_listings array
if (.not. associated(ace_tables_)) then
message = "No ACE table listings present in cross_sections.xml file!"
call fatal_error()
message = "No ACE table listings present in cross_sections.xml file!"
call fatal_error()
else
n_listings = size(ace_tables_)
allocate(xs_listings(n_listings))
n_listings = size(ace_tables_)
allocate(xs_listings(n_listings))
end if
do i = 1, n_listings
listing => xs_listings(i)
listing => xs_listings(i)
! copy a number of attributes
listing % name = trim(ace_tables_(i) % name)
listing % alias = trim(ace_tables_(i) % alias)
listing % zaid = ace_tables_(i) % zaid
listing % awr = ace_tables_(i) % awr
listing % kT = ace_tables_(i) % temperature
listing % location = ace_tables_(i) % location
! copy a number of attributes
listing % name = trim(ace_tables_(i) % name)
listing % alias = trim(ace_tables_(i) % alias)
listing % zaid = ace_tables_(i) % zaid
listing % awr = ace_tables_(i) % awr
listing % kT = ace_tables_(i) % temperature
listing % location = ace_tables_(i) % location
! determine type of cross section
if (ends_with(listing % name, 'c')) then
! determine type of cross section
if (ends_with(listing % name, 'c')) then
listing % type = ACE_NEUTRON
elseif (ends_with(listing % name, 't')) then
elseif (ends_with(listing % name, 't')) then
listing % type = ACE_THERMAL
end if
end if
! set filetype, record length, and number of entries
listing % filetype = filetype
listing % recl = recl
listing % entries = entries
! set filetype, record length, and number of entries
listing % filetype = filetype
listing % recl = recl
listing % entries = entries
! determine metastable state
if (ace_tables_(i) % metastable == 0) then
! determine metastable state
if (ace_tables_(i) % metastable == 0) then
listing % metastable = .false.
else
else
listing % metastable = .true.
end if
end if
! determine path of cross section table
if (starts_with(ace_tables_(i) % path, '/')) then
! determine path of cross section table
if (starts_with(ace_tables_(i) % path, '/')) then
listing % path = ace_tables_(i) % path
else
else
if (ends_with(directory,'/')) then
listing % path = trim(directory) // trim(ace_tables_(i) % path)
listing % path = trim(directory) // trim(ace_tables_(i) % path)
else
listing % path = trim(directory) // '/' // trim(ace_tables_(i) % path)
listing % path = trim(directory) // '/' // trim(ace_tables_(i) % path)
end if
end if
end if
! create dictionary entry for both name and alias
call dict_add_key(xs_listing_dict, listing % name, i)
call dict_add_key(xs_listing_dict, listing % alias, i)
! create dictionary entry for both name and alias
call dict_add_key(xs_listing_dict, listing % name, i)
call dict_add_key(xs_listing_dict, listing % alias, i)
end do
end subroutine read_cross_sections_xml

6
src/main.F90
View file

@ -1,7 +1,7 @@
program main
use constants
use criticality, only: run_criticality
use eigenvalue, only: run_eigenvalue
use finalize, only: finalize_run
use fixed_source, only: run_fixedsource
use global
@ -17,8 +17,8 @@ program main
select case (run_mode)
case (MODE_FIXEDSOURCE)
call run_fixedsource()
case (MODE_CRITICALITY)
call run_criticality()
case (MODE_EIGENVALUE)
call run_eigenvalue()
case (MODE_PLOTTING)
call run_plot()
case (MODE_TALLIES)

4
src/mesh.F90
View file

@ -93,7 +93,7 @@ contains
!===============================================================================
! MESH_INDICES_TO_BIN maps (i,j) or (i,j,k) indices to a single bin number for
! use in a TallyObject scores array
! use in a TallyObject results array
!===============================================================================
function mesh_indices_to_bin(m, ijk, surface_current) result(bin)
@ -126,7 +126,7 @@ contains
end function mesh_indices_to_bin
!===============================================================================
! BIN_TO_MESH_INDICES maps a single mesh bin from a TallyObject scores array to
! BIN_TO_MESH_INDICES maps a single mesh bin from a TallyObject results array to
! (i,j) or (i,j,k) indices
!===============================================================================

271
src/output.F90
View file

@ -32,35 +32,35 @@ contains
subroutine title()
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'
' .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(UNIT=OUTPUT_UNIT, FMT=*) &
' Developed At: Massachusetts Institute of Technology'
' Developed At: Massachusetts Institute of Technology'
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",7X,I1,".",I1,".",I1)') &
VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
#ifdef GIT_SHA1
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",6X,A)') GIT_SHA1
#endif
! Write the date and time
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A)') &
time_stamp()
time_stamp()
#ifdef MPI
! Write number of processors
write(UNIT=OUTPUT_UNIT, FMT='(6X,"MPI Processes:",1X,A)') &
trim(to_str(n_procs))
trim(to_str(n_procs))
#endif
end subroutine title
@ -77,7 +77,7 @@ contains
call date_and_time(DATE=date_, TIME=time_)
current_time = date_(1:4) // "-" // date_(5:6) // "-" // date_(7:8) // &
" " // time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6)
" " // time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6)
end function time_stamp
@ -125,13 +125,13 @@ contains
select case (header_level)
case (1)
write(UNIT=unit_, FMT='(/3(1X,A/))') repeat('=', 75), &
repeat('=', n) // '> ' // trim(line) // ' <' // &
repeat('=', m), repeat('=', 75)
repeat('=', n) // '> ' // trim(line) // ' <' // &
repeat('=', m), repeat('=', 75)
case (2)
write(UNIT=unit_, FMT='(/2(1X,A/))') trim(line), repeat('-', 75)
case (3)
write(UNIT=unit_, FMT='(/1X,A/)') repeat('=', n) // '> ' // &
trim(line) // ' <' // repeat('=', m)
trim(line) // ' <' // repeat('=', m)
end select
end subroutine header
@ -145,11 +145,11 @@ contains
if (master) then
write(UNIT=OUTPUT_UNIT, FMT='(1X,A,1X,I1,".",I1,".",I1)') &
"OpenMC version", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
"OpenMC version", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
write(UNIT=OUTPUT_UNIT, FMT=*) "Copyright (c) 2011-2012 &
&Massachusetts Institute of Technology"
&Massachusetts Institute of Technology"
write(UNIT=OUTPUT_UNIT, FMT=*) "MIT/X license at &
&<http://mit-crpg.github.com/openmc/license.html>"
&<http://mit-crpg.github.com/openmc/license.html>"
end if
end subroutine print_version
@ -207,7 +207,7 @@ contains
else
! Determine last space in current line
i_end = i_start + index(message(i_start+1:i_start+line_wrap), &
' ', BACK=.true.)
' ', BACK=.true.)
! Write up to last space
write(ou, fmt='(1X,A)') message(i_start+1:i_end-1)
@ -270,8 +270,8 @@ contains
l => lattices(coord % lattice)
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)) // ')'
p % coord % lattice_x)) // ',' // trim(to_str(&
p % coord % lattice_y)) // ')'
end if
! Print local coordinates
@ -390,7 +390,7 @@ contains
case default
index_surf = abs(c % surfaces(i))
string = trim(string) // ' ' // to_str(sign(&
surfaces(index_surf) % id, c % surfaces(i)))
surfaces(index_surf) % id, c % surfaces(i)))
end select
end do
write(unit_,*) ' Surface Specification:' // trim(string)
@ -484,6 +484,7 @@ contains
integer :: i ! loop index for coefficients
integer :: unit_ ! unit to write to
character(MAX_LINE_LEN) :: string
type(Cell), pointer :: c => null()
! set default unit if not specified
if (present(unit)) then
@ -539,7 +540,9 @@ contains
string = ""
if (allocated(surf % neighbor_pos)) then
do i = 1, size(surf % neighbor_pos)
string = trim(string) // ' ' // to_str(surf % neighbor_pos(i))
c => cells(abs(surf % neighbor_pos(i)))
string = trim(string) // ' ' // to_str(&
sign(c % id, surf % neighbor_pos(i)))
end do
end if
write(unit_,*) ' Positive Neighbors = ' // trim(string)
@ -548,7 +551,9 @@ contains
string = ""
if (allocated(surf % neighbor_neg)) then
do i = 1, size(surf % neighbor_neg)
string = trim(string) // ' ' // to_str(surf % neighbor_neg(i))
c => cells(abs(surf % neighbor_neg(i)))
string = trim(string) // ' ' // to_str(&
sign(c % id, surf % neighbor_neg(i)))
end do
end if
write(unit_,*) ' Negative Neighbors =' // trim(string)
@ -595,7 +600,7 @@ contains
! Write total atom density in atom/b-cm
write(unit_,*) ' Atom Density = ' // trim(to_str(mat % density)) &
// ' atom/b-cm'
// ' atom/b-cm'
! Write atom density for each nuclide in material
write(unit_,*) ' Nuclides:'
@ -603,7 +608,7 @@ contains
nuc => nuclides(mat % nuclide(i))
density = mat % atom_density(i)
string = ' ' // trim(nuc % name) // ' = ' // &
trim(to_str(density)) // ' atom/b-cm'
trim(to_str(density)) // ' atom/b-cm'
write(unit_,*) trim(string)
end do
@ -742,7 +747,7 @@ contains
string = ""
do i = 1, t % filters(j) % n_bins + 1
string = trim(string) // ' ' // trim(to_str(&
t % filters(j) % real_bins(i)))
t % filters(j) % real_bins(i)))
end do
write(unit_,*) ' Incoming Energy Bins:' // trim(string)
end if
@ -753,7 +758,7 @@ contains
string = ""
do i = 1, t % filters(j) % n_bins + 1
string = trim(string) // ' ' // trim(to_str(&
t % filters(j) % real_bins(i)))
t % filters(j) % real_bins(i)))
end do
write(unit_,*) ' Outgoing Energy Bins:' // trim(string)
end if
@ -765,10 +770,10 @@ contains
write(unit_,fmt='(A)',advance='no') ' total'
else
write(unit_,fmt='(A)',advance='no') ' ' // trim(adjustl(&
nuclides(t % nuclide_bins(i)) % name))
nuclides(t % nuclide_bins(i)) % name))
end if
if (mod(i,4) == 0 .and. i /= t % n_nuclide_bins) &
write(unit_,'(/18X)',advance='no')
write(unit_,'(/18X)',advance='no')
end do
write(unit_,*)
@ -938,8 +943,8 @@ contains
end if
write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,A4,1X,I6,1X,I11,1X,I11)') &
reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, &
law(1:4), rxn % threshold, size_angle, size_energy
reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, &
law(1:4), rxn % threshold, size_angle, size_energy
! Accumulate data size
size_xs = size_xs + (nuc % n_grid - rxn%threshold + 1) * 8
@ -976,11 +981,11 @@ contains
write(unit_,*) ' Memory Requirements'
write(unit_,*) ' Cross sections = ' // trim(to_str(size_xs)) // ' bytes'
write(unit_,*) ' Secondary angle distributions = ' // &
trim(to_str(size_angle_total)) // ' bytes'
trim(to_str(size_angle_total)) // ' bytes'
write(unit_,*) ' Secondary energy distributions = ' // &
trim(to_str(size_energy_total)) // ' bytes'
trim(to_str(size_energy_total)) // ' bytes'
write(unit_,*) ' Probability Tables = ' // &
trim(to_str(size_urr)) // ' bytes'
trim(to_str(size_urr)) // ' bytes'
write(unit_,*) ' Total = ' // trim(to_str(size_total)) // ' bytes'
! Blank line at end of nuclide
@ -1016,28 +1021,28 @@ contains
! Inelastic data
write(unit_,*) ' # of Incoming Energies (Inelastic) = ' // &
trim(to_str(sab % n_inelastic_e_in))
trim(to_str(sab % n_inelastic_e_in))
write(unit_,*) ' # of Outgoing Energies (Inelastic) = ' // &
trim(to_str(sab % n_inelastic_e_out))
trim(to_str(sab % n_inelastic_e_out))
write(unit_,*) ' # of Outgoing Angles (Inelastic) = ' // &
trim(to_str(sab % n_inelastic_mu))
trim(to_str(sab % n_inelastic_mu))
write(unit_,*) ' Threshold for Inelastic = ' // &
trim(to_str(sab % threshold_inelastic))
trim(to_str(sab % threshold_inelastic))
! Elastic data
if (sab % n_elastic_e_in > 0) then
write(unit_,*) ' # of Incoming Energies (Elastic) = ' // &
trim(to_str(sab % n_elastic_e_in))
trim(to_str(sab % n_elastic_e_in))
write(unit_,*) ' # of Outgoing Angles (Elastic) = ' // &
trim(to_str(sab % n_elastic_mu))
trim(to_str(sab % n_elastic_mu))
write(unit_,*) ' Threshold for Elastic = ' // &
trim(to_str(sab % threshold_elastic))
trim(to_str(sab % threshold_elastic))
end if
! Determine memory used by S(a,b) table and write out
size_sab = 8 * (sab % n_inelastic_e_in * (2 + sab % n_inelastic_e_out * &
(1 + sab % n_inelastic_mu)) + sab % n_elastic_e_in * &
(2 + sab % n_elastic_mu))
(1 + sab % n_inelastic_mu)) + sab % n_elastic_e_in * &
(2 + sab % n_elastic_mu))
write(unit_,*) ' Memory Used = ' // trim(to_str(size_sab)) // ' bytes'
! Blank line at end
@ -1065,31 +1070,31 @@ contains
call header("OpenMC Monte Carlo Code", unit=UNIT_SUMMARY, level=1)
write(UNIT=UNIT_SUMMARY, FMT=*) &
"Copyright: 2011-2012 Massachusetts Institute of Technology"
"Copyright: 2011-2012 Massachusetts Institute of Technology"
write(UNIT=UNIT_SUMMARY, FMT='(1X,A,7X,2(I1,"."),I1)') &
"Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
"Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
#ifdef GIT_SHA1
write(UNIT=UNIT_SUMMARY, FMT='(1X,"Git SHA1:",6X,A)') GIT_SHA1
#endif
write(UNIT=UNIT_SUMMARY, FMT='(1X,"Date/Time:",5X,A)') &
time_stamp()
time_stamp()
! Write information on number of processors
#ifdef MPI
write(UNIT=UNIT_SUMMARY, FMT='(1X,"MPI Processes:",1X,A)') &
trim(to_str(n_procs))
trim(to_str(n_procs))
#endif
! Display problem summary
call header("PROBLEM SUMMARY", unit=UNIT_SUMMARY)
select case(run_mode)
case (MODE_CRITICALITY)
write(UNIT_SUMMARY,100) 'Problem type:', 'Criticality'
case (MODE_EIGENVALUE)
write(UNIT_SUMMARY,100) 'Problem type:', 'k eigenvalue'
write(UNIT_SUMMARY,101) 'Number of Batches:', n_batches
write(UNIT_SUMMARY,101) 'Number of Inactive Batches:', n_inactive
write(UNIT_SUMMARY,101) 'Generations per Batch:', gen_per_batch
case (MODE_FIXEDSOURCE)
write(UNIT_SUMMARY,100) 'Problem type:', 'External Source'
write(UNIT_SUMMARY,100) 'Problem type:', 'fixed source'
end select
write(UNIT_SUMMARY,101) 'Number of Particles:', n_particles
@ -1122,7 +1127,7 @@ contains
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"
n_grid*n_nuclides_total*4)) // " bytes"
! print summary of variance reduction
call header("VARIANCE REDUCTION", unit=UNIT_SUMMARY)
@ -1232,16 +1237,16 @@ contains
! write out information batch and option independent output
write(UNIT=OUTPUT_UNIT, FMT='(2X,I5)', ADVANCE='NO') current_batch
write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') &
k_batch(current_batch)
k_batch(current_batch)
! write out entropy info
if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
entropy(current_batch)
entropy(current_batch)
! write out accumulated k-effective if after first active batch
if (current_batch > n_inactive + 1) then
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') &
keff, keff_std
keff, keff_std
else
write(UNIT=OUTPUT_UNIT, FMT='(23X)', ADVANCE='NO')
end if
@ -1249,11 +1254,11 @@ contains
! write out cmfd keff if it is active
if (cmfd_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % keff
cmfd % keff
! write out cmfd entopy
if (cmfd_on .and. entropy_on) write(UNIT=OUTPUT_UNIT, &
FMT='(3X, F8.5)', ADVANCE='NO') cmfd % entropy
FMT='(3X, F8.5)', ADVANCE='NO') cmfd % entropy
! next line
write(UNIT=OUTPUT_UNIT, FMT=*)
@ -1280,18 +1285,18 @@ contains
! Write plotting origin
write(ou,100) "Origin:", trim(to_str(pl % origin(1))) // &
" " // trim(to_str(pl % origin(2))) // " " // &
trim(to_str(pl % origin(3)))
" " // trim(to_str(pl % origin(2))) // " " // &
trim(to_str(pl % origin(3)))
! Write plotting width
if (pl % type == PLOT_TYPE_SLICE) then
write(ou,100) "Width:", trim(to_str(pl % width(1))) // &
" " // trim(to_str(pl % width(2)))
" " // trim(to_str(pl % width(2)))
write(ou,100) "Coloring:", trim(to_str(pl % color_by))
write(ou,100) "Basis:", trim(to_str(pl % basis))
write(ou,100) "Pixels:", trim(to_str(pl % pixels(1))) // " " // &
trim(to_str(pl % pixels(2)))
trim(to_str(pl % pixels(2)))
end if
write(ou,*)
@ -1324,7 +1329,7 @@ contains
write(ou,100) " Reading cross sections", time_read_xs % elapsed
write(ou,100) " Unionizing energy grid", time_unionize % elapsed
write(ou,100) "Total time in simulation", time_inactive % elapsed + &
time_active % elapsed
time_active % elapsed
write(ou,100) " Time in transport only", time_transport % elapsed
if(cmfd_run) write(ou,100) "Total CMFD time", time_cmfd % elapsed
write(ou,100) " Time in inactive batches", time_inactive % elapsed
@ -1338,14 +1343,14 @@ contains
if (restart_run) then
total_particles = n_particles * (n_batches - &
restart_batch) * gen_per_batch
restart_batch) * gen_per_batch
else
total_particles = n_particles * n_batches * gen_per_batch
end if
! display calculate rate
speed = real(total_particles) / (time_inactive % elapsed + &
time_active % elapsed)
time_active % elapsed)
string = to_str(speed)
write(ou,101) "Calculation Rate", trim(string)
@ -1363,13 +1368,13 @@ contains
! write global tallies
write(ou,102) "k-effective (Analog)", global_tallies(K_ANALOG) % sum, &
global_tallies(K_ANALOG) % sum_sq
global_tallies(K_ANALOG) % sum_sq
write(ou,102) "k-effective (Collision)", global_tallies(K_COLLISION) % sum, &
global_tallies(K_COLLISION) % sum_sq
global_tallies(K_COLLISION) % sum_sq
write(ou,102) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) % sum, &
global_tallies(K_TRACKLENGTH) % sum_sq
global_tallies(K_TRACKLENGTH) % sum_sq
write(ou,102) "Leakage Fraction", global_tallies(LEAKAGE) % sum, &
global_tallies(LEAKAGE) % sum_sq
global_tallies(LEAKAGE) % sum_sq
write(ou,*)
! format for write statements
@ -1393,11 +1398,11 @@ contains
integer :: l ! loop index for nonPN scores
integer :: type ! type of tally filter
integer :: indent ! number of spaces to preceed output
integer :: filter_index ! index in scores array for filters
integer :: filter_index ! index in results array for filters
integer :: score_index ! scoring bin index
integer :: i_nuclide ! index in nuclides array
integer :: i_listing ! index in xs_listings array
integer :: n_order ! loop index for scattering orders
integer :: n_order ! loop index for scattering orders
real(8) :: t_value ! t-values for confidence intervals
real(8) :: alpha ! significance level for CI
character(MAX_FILE_LEN) :: filename ! name of output file
@ -1459,14 +1464,14 @@ contains
! Multiply uncertainty by t-value
if (confidence_intervals) then
do k = 1, size(t % scores, 2)
do j = 1, size(t % scores, 1)
do k = 1, size(t % results, 2)
do j = 1, size(t % results, 1)
! Calculate t-value for confidence intervals
if (confidence_intervals) then
alpha = ONE - CONFIDENCE_LEVEL
t_value = t_percentile(ONE - alpha/TWO, t % n_realizations - 1)
end if
t % scores(j,k) % sum_sq = t_value * t % scores(j,k) % sum_sq
t % results(j,k) % sum_sq = t_value * t % results(j,k) % sum_sq
end do
end do
end if
@ -1474,10 +1479,10 @@ contains
! Write header block
if (t % label == "") then
call header("TALLY " // trim(to_str(t % id)), unit=UNIT_TALLY, &
level=3)
level=3)
else
call header("TALLY " // trim(to_str(t % id)) // ": " &
// trim(t % label), unit=UNIT_TALLY, level=3)
// trim(t % label), unit=UNIT_TALLY, level=3)
endif
! Handle surface current tallies separately
@ -1486,7 +1491,7 @@ contains
cycle
end if
! WARNING: Admittedly, the logic for moving for printing scores is
! WARNING: Admittedly, the logic for moving for printing results is
! extremely confusing and took quite a bit of time to get correct. The
! logic is structured this way since it is not practical to have a do
! loop for each filter variable (given that only a few filters are likely
@ -1517,7 +1522,7 @@ contains
indent = indent - 2
! =================================================================
! VALID BIN -- WRITE FILTER INFORMATION OR EXIT TO WRITE SCORES
! VALID BIN -- WRITE FILTER INFORMATION OR EXIT TO WRITE RESULTS
else
! Check if this is last filter
@ -1526,7 +1531,7 @@ contains
! Print current filter information
type = t % filters(j) % type
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(filter_name(type)), trim(get_label(t, j))
trim(filter_name(type)), trim(get_label(t, j))
indent = indent + 2
j = j + 1
end if
@ -1537,7 +1542,7 @@ contains
if (t % n_filters > 0) then
type = t % filters(j) % type
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(filter_name(type)), trim(get_label(t, j))
trim(filter_name(type)), trim(get_label(t, j))
end if
! Determine scoring index for this bin combination -- note that unlike
@ -1549,8 +1554,8 @@ contains
else
filter_index = 1
end if
! Write scores for this filter bin combination
! Write results for this filter bin combination
score_index = 0
if (t % n_filters > 0) indent = indent + 2
do n = 1, t % n_nuclide_bins
@ -1558,11 +1563,11 @@ contains
i_nuclide = t % nuclide_bins(n)
if (i_nuclide == -1) then
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
"Total Material"
"Total Material"
else
i_listing = nuclides(i_nuclide) % listing
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(xs_listings(i_listing) % alias)
trim(xs_listings(i_listing) % alias)
end if
indent = indent + 2
@ -1639,7 +1644,7 @@ contains
integer :: n ! number of incoming energy bins
integer :: len1 ! length of string
integer :: len2 ! length of string
integer :: filter_index ! index in scores array for filters
integer :: filter_index ! index in results array for filters
logical :: print_ebin ! should incoming energy bin be displayed?
character(MAX_LINE_LEN) :: string
type(StructuredMesh), pointer :: m => null()
@ -1680,110 +1685,110 @@ contains
! Write incoming energy bin
write(UNIT=UNIT_TALLY, FMT='(3X,A,1X,A)') &
"Incoming Energy", trim(get_label(t, i_filter_ein))
"Incoming Energy", trim(get_label(t, i_filter_ein))
end if
! Left Surface
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.)
mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Left", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Outgoing Current to Left", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Left", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Incoming Current from Left", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
! Right Surface
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Right", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Incoming Current from Right", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Right", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Outgoing Current to Right", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
! Back Surface
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.)
mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Back", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Outgoing Current to Back", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_FRONT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Back", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Incoming Current from Back", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
! Front Surface
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Front", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Incoming Current from Front", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_FRONT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Front", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Outgoing Current to Front", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
! Bottom Surface
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.)
mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Bottom", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Outgoing Current to Bottom", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_TOP
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Bottom", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Incoming Current from Bottom", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
! Top Surface
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Top", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Incoming Current from Top", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_TOP
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Top", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Outgoing Current to Top", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
end do
end do
@ -1831,10 +1836,10 @@ contains
call bin_to_mesh_indices(m, bin, ijk)
if (m % n_dimension == 2) then
label = "Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ")"
trim(to_str(ijk(2))) // ")"
elseif (m % n_dimension == 3) then
label = "Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
end if
case (FILTER_ENERGYIN, FILTER_ENERGYOUT)
E0 = t % filters(i_filter) % real_bins(bin)

22
src/physics.F90
View file

@ -19,7 +19,7 @@ module physics
use search, only: binary_search
use string, only: to_str
use tally, only: score_analog_tally, score_tracklength_tally, &
score_surface_current, add_to_score
score_surface_current
implicit none
@ -103,12 +103,10 @@ contains
if (associated(active_tracklength_tallies)) &
call score_tracklength_tally(distance)
! Score global tracklength tallies
if (tallies_on) then
! Score track-length estimate of k-eff
call add_to_score(global_tallies(K_TRACKLENGTH), &
p % wgt * distance * material_xs % nu_fission)
end if
! Score track-length estimate of k-eff
if (tallies_on) global_tallies(K_TRACKLENGTH) % value = &
global_tallies(K_TRACKLENGTH) % value + p % wgt * distance * &
material_xs % nu_fission
if (d_collision > d_boundary) then
! ====================================================================
@ -132,10 +130,9 @@ contains
! PARTICLE HAS COLLISION
! Score collision estimate of keff
if (tallies_on) then
call add_to_score(global_tallies(K_COLLISION), &
p % wgt * material_xs % nu_fission / material_xs % total)
end if
if (tallies_on) global_tallies(K_COLLISION) % value = &
global_tallies(K_COLLISION) % value + p % wgt * &
material_xs % nu_fission / material_xs % total
p % surface = NONE
call collision()
@ -908,7 +905,8 @@ contains
end if
! Add to analog estimate of keff
call add_to_score(global_tallies(K_ANALOG), nu/weight * keff)
global_tallies(K_ANALOG) % value = &
global_tallies(K_ANALOG) % value + nu/weight * keff
! Bank source neutrons
if (nu == 0 .or. n_bank == 3*work) return

2
src/plot.F90
View file

@ -95,7 +95,7 @@ contains
allocate(p)
call initialize_particle()
p % coord % xyz = xyz
p % coord % uvw = (/ 1, 0, 0 /)
p % coord % uvw = (/ 0.5, 0.5, 0.5 /)
p % coord % universe = BASE_UNIVERSE
do y = 1, img % height

2
src/source_header.F90
View file

@ -4,7 +4,7 @@ module source_header
!===============================================================================
! EXTSOURCE describes an external source of neutrons for a fixed-source problem
! or for the starting source in a criticality problem
! or for the starting source in a k eigenvalue problem
!===============================================================================
type ExtSource

124
src/state_point.F90
View file

@ -82,7 +82,7 @@ contains
! If using the no-tally-reduction method, we need to collect tally
! results before writing them to the state point file.
call write_tally_scores_nr(fh)
call write_tally_results_nr(fh)
elseif (master) then
! Write number of realizations
@ -93,7 +93,7 @@ contains
call MPI_FILE_WRITE(fh, N_GLOBAL_TALLIES, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpi_err)
call MPI_FILE_WRITE(fh, global_tallies, N_GLOBAL_TALLIES, &
MPI_TALLYSCORE, MPI_STATUS_IGNORE, mpi_err)
MPI_TALLYRESULT, MPI_STATUS_IGNORE, mpi_err)
if (tallies_on) then
! Indicate that tallies are on
@ -101,14 +101,14 @@ contains
call MPI_FILE_WRITE(fh, temp, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpi_err)
! Write all tally scores
TALLY_SCORES: do i = 1, n_tallies
! Write all tally results
TALLY_RESULTS: do i = 1, n_tallies
t => tallies(i)
n = size(t % scores, 1) * size(t % scores, 2)
call MPI_FILE_WRITE(fh, t % scores, n, MPI_TALLYSCORE, &
n = size(t % results, 1) * size(t % results, 2)
call MPI_FILE_WRITE(fh, t % results, n, MPI_TALLYRESULT, &
MPI_STATUS_IGNORE, mpi_err)
end do TALLY_SCORES
end do TALLY_RESULTS
else
! Indicate that tallies are off
temp = 0
@ -120,7 +120,7 @@ contains
! ==========================================================================
! SOURCE BANK
if (run_mode == MODE_CRITICALITY) then
if (run_mode == MODE_EIGENVALUE) then
if (source_separate) then
! If the user has specified that the source sites should be written in
! a separate file, we make a call to the appropriate subroutine to
@ -179,8 +179,8 @@ contains
! Write out current batch number
write(UNIT_STATE) current_batch
! Write out information for criticality run
if (run_mode == MODE_CRITICALITY) then
! Write out information for eigenvalue run
if (run_mode == MODE_EIGENVALUE) then
write(UNIT_STATE) n_inactive, gen_per_batch
write(UNIT_STATE) k_batch(1:current_batch)
write(UNIT_STATE) entropy(1:current_batch)
@ -209,7 +209,7 @@ contains
! Number of realizations
write(UNIT_STATE) t % n_realizations
! Write size of each dimension of tally scores array
! Write size of each dimension of tally results array
write(UNIT_STATE) t % total_score_bins
write(UNIT_STATE) t % total_filter_bins
@ -263,25 +263,25 @@ contains
! Indicate that tallies are on
write(UNIT_STATE) 1
TALLY_SCORES: do i = 1, n_tallies
TALLY_RESULTS: do i = 1, n_tallies
! Get pointer to tally
t => tallies(i)
! Write tally sum and sum_sq for each bin
do k = 1, size(t % scores, 2)
do j = 1, size(t % scores, 1)
write(UNIT_STATE) t % scores(j,k) % sum
write(UNIT_STATE) t % scores(j,k) % sum_sq
do k = 1, size(t % results, 2)
do j = 1, size(t % results, 1)
write(UNIT_STATE) t % results(j,k) % sum
write(UNIT_STATE) t % results(j,k) % sum_sq
end do
end do
end do TALLY_SCORES
end do TALLY_RESULTS
else
! Indicate that tallies are off
write(UNIT_STATE) 0
end if
! Write out source bank
if (run_mode == MODE_CRITICALITY) then
if (run_mode == MODE_EIGENVALUE) then
if (source_separate) then
! If the user has specified that the source sites should be written in
! a separate file, we make a call to the appropriate subroutine to
@ -349,8 +349,8 @@ contains
call MPI_FILE_WRITE(fh, current_batch, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpi_err)
! Write out information for criticality run
if (run_mode == MODE_CRITICALITY) then
! Write out information for eigenvalue run
if (run_mode == MODE_EIGENVALUE) then
call MPI_FILE_WRITE(fh, n_inactive, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpi_err)
call MPI_FILE_WRITE(fh, gen_per_batch, 1, MPI_INTEGER, &
@ -454,10 +454,10 @@ contains
#ifdef MPI
!===============================================================================
! WRITE_TALLY_SCORES_NR
! WRITE_TALLY_RESULTS_NR
!===============================================================================
subroutine write_tally_scores_nr(fh)
subroutine write_tally_results_nr(fh)
integer, intent(in) :: fh ! file handle
@ -466,7 +466,7 @@ contains
integer :: m ! number of score bins
integer :: temp ! temporary variable
integer :: n_bins ! total number of bins
real(8), allocatable :: tally_temp(:,:,:) ! contiguous array of scores
real(8), allocatable :: tally_temp(:,:,:) ! contiguous array of results
real(8) :: global_temp(2,N_GLOBAL_TALLIES)
real(8) :: dummy ! temporary receive buffer for non-root reduces
type(TallyObject), pointer :: t => null()
@ -518,20 +518,20 @@ contains
MPI_STATUS_IGNORE, mpi_err)
end if
! Write all tally scores
TALLY_SCORES: do i = 1, n_tallies
! Write all tally results
TALLY_RESULTS: do i = 1, n_tallies
t => tallies(i)
! Determine size of tally scores array
m = size(t % scores, 1)
n = size(t % scores, 2)
! Determine size of tally results array
m = size(t % results, 1)
n = size(t % results, 2)
n_bins = m*n*2
! Allocate array for storing sums and sums of squares, but
! contiguously in memory for each
allocate(tally_temp(2,m,n))
tally_temp(1,:,:) = t % scores(:,:) % sum
tally_temp(2,:,:) = t % scores(:,:) % sum_sq
tally_temp(1,:,:) = t % results(:,:) % sum
tally_temp(2,:,:) = t % results(:,:) % sum_sq
if (master) then
! The MPI_IN_PLACE specifier allows the master to copy values into
@ -544,10 +544,10 @@ contains
MPI_STATUS_IGNORE, mpi_err)
! At the end of the simulation, store the results back in the
! regular TallyScores array
! regular TallyResults array
if (current_batch == n_batches) then
t % scores(:,:) % sum = tally_temp(1,:,:)
t % scores(:,:) % sum_sq = tally_temp(2,:,:)
t % results(:,:) % sum = tally_temp(1,:,:)
t % results(:,:) % sum_sq = tally_temp(2,:,:)
end if
else
! Receive buffer not significant at other processors
@ -557,7 +557,7 @@ contains
! Deallocate temporary copy of tally results
deallocate(tally_temp)
end do TALLY_SCORES
end do TALLY_RESULTS
else
if (master) then
! Indicate that tallies are off
@ -567,7 +567,7 @@ contains
end if
end if
end subroutine write_tally_scores_nr
end subroutine write_tally_results_nr
#endif
!===============================================================================
@ -644,8 +644,8 @@ contains
call MPI_FILE_READ_ALL(fh, restart_batch, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpi_err)
! Read information specific to criticality run
if (mode == MODE_CRITICALITY) then
! Read information specific to eigenvalue run
if (mode == MODE_EIGENVALUE) then
call MPI_FILE_READ_ALL(fh, n_inactive, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpi_err)
call MPI_FILE_READ_ALL(fh, gen_per_batch, 1, MPI_INTEGER, &
@ -687,12 +687,12 @@ contains
call MPI_FILE_READ(fh, tallies(i) % n_realizations, 1, &
MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err)
! Read dimensions of tally filters and scores and make sure they
! Read dimensions of tally filters and results and make sure they
! match
call MPI_FILE_READ(fh, temp, 2, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpi_err)
if (temp(1) /= size(tallies(i) % scores, 1) .or. &
temp(2) /= size(tallies(i) % scores, 2)) then
if (temp(1) /= size(tallies(i) % results, 1) .or. &
temp(2) /= size(tallies(i) % results, 2)) then
message = "Tally dimensions do not match in state point."
call fatal_error()
end if
@ -736,7 +736,7 @@ contains
MPI_STATUS_IGNORE, mpi_err)
deallocate(int_array)
! Read number of scores
! Read number of results
call MPI_FILE_READ(fh, temp, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpi_err)
@ -760,7 +760,7 @@ contains
! Read global tally data
call MPI_FILE_READ(fh, global_tallies, N_GLOBAL_TALLIES, &
MPI_TALLYSCORE, MPI_STATUS_IGNORE, mpi_err)
MPI_TALLYRESULT, MPI_STATUS_IGNORE, mpi_err)
! Check if tally results are present
call MPI_FILE_READ(fh, temp, 1, MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err)
@ -770,18 +770,18 @@ contains
! Read sum and sum squared
if (temp(1) == 1) then
TALLY_SCORES: do i = 1, n_tallies
n = size(tallies(i) % scores, 1) * size(tallies(i) % scores, 2)
call MPI_FILE_READ(fh, tallies(i) % scores, n, MPI_TALLYSCORE, &
TALLY_RESULTS: do i = 1, n_tallies
n = size(tallies(i) % results, 1) * size(tallies(i) % results, 2)
call MPI_FILE_READ(fh, tallies(i) % results, n, MPI_TALLYRESULT, &
MPI_STATUS_IGNORE, mpi_err)
end do TALLY_SCORES
end do TALLY_RESULTS
end if
end if
! ==========================================================================
! SOURCE BANK
if (run_mode == MODE_CRITICALITY) then
if (run_mode == MODE_EIGENVALUE) then
! Get current offset for master
if (master) call MPI_FILE_GET_POSITION(fh, offset, mpi_err)
@ -841,8 +841,8 @@ contains
! Read batch number to restart at
read(UNIT_STATE) restart_batch
! Read information specific to criticality run
if (mode == MODE_CRITICALITY) then
! Read information specific to eigenvalue run
if (mode == MODE_EIGENVALUE) then
read(UNIT_STATE) n_inactive, gen_per_batch
read(UNIT_STATE) k_batch(1:restart_batch)
read(UNIT_STATE) entropy(1:restart_batch)
@ -886,11 +886,11 @@ contains
! Read number of realizations
read(UNIT_STATE) tallies(i) % n_realizations
! Read dimensions of tally filters and scores and make sure they
! Read dimensions of tally filters and results and make sure they
! match
read(UNIT_STATE) temp(1:2)
if (temp(1) /= size(tallies(i) % scores, 1) .or. &
temp(2) /= size(tallies(i) % scores, 2)) then
if (temp(1) /= size(tallies(i) % results, 1) .or. &
temp(2) /= size(tallies(i) % results, 2)) then
message = "Tally dimensions do not match in state point."
call fatal_error()
end if
@ -927,7 +927,7 @@ contains
read(UNIT_STATE) int_array
deallocate(int_array)
! Read number of scores
! Read number of results
read(UNIT_STATE) temp(1)
! Read nuclide bins
@ -955,19 +955,19 @@ contains
! Read sum and sum squared
read(UNIT_STATE) temp(1)
if (temp(1) == 1) then
TALLY_SCORES: do i = 1, n_tallies
do k = 1, size(tallies(i) % scores, 2)
do j = 1, size(tallies(i) % scores, 1)
read(UNIT_STATE) tallies(i) % scores(j,k) % sum
read(UNIT_STATE) tallies(i) % scores(j,k) % sum_sq
TALLY_RESULTS: do i = 1, n_tallies
do k = 1, size(tallies(i) % results, 2)
do j = 1, size(tallies(i) % results, 1)
read(UNIT_STATE) tallies(i) % results(j,k) % sum
read(UNIT_STATE) tallies(i) % results(j,k) % sum_sq
end do
end do
end do TALLY_SCORES
end do TALLY_RESULTS
end if
end if
! Read source bank for criticality run
if (mode == MODE_CRITICALITY .and. run_mode /= MODE_TALLIES) then
! Read source bank for eigenvalue run
if (mode == MODE_EIGENVALUE .and. run_mode /= MODE_TALLIES) then
read(UNIT_STATE) source_bank
end if
@ -995,7 +995,7 @@ contains
call write_message(1)
end if
! For criticality calculations, turn on tallies if we've reached active
! For eigenvalue calculations, turn on tallies if we've reached active
! batches
if (current_batch == n_inactive) then
tallies_on = .true.

272
src/tally.F90
View file

@ -6,12 +6,12 @@ module tally
use global
use math, only: t_percentile, calc_pn
use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices, &
mesh_indices_to_bin, mesh_intersects
mesh_indices_to_bin, mesh_intersects
use mesh_header, only: StructuredMesh
use output, only: header
use search, only: binary_search
use string, only: to_str
use tally_header, only: TallyScore, TallyMapItem, TallyMapElement
use tally_header, only: TallyResult, TallyMapItem, TallyMapElement
#ifdef MPI
use mpi
@ -76,10 +76,10 @@ contains
end if
! =======================================================================
! CALCULATE SCORES AND ACCUMULATE TALLY
! CALCULATE RESULTS AND ACCUMULATE TALLY
! If we have made it here, we have a scoring combination of bins for this
! tally -- now we need to determine where in the scores array we should
! tally -- now we need to determine where in the results array we should
! be accumulating the tally values
! Determine scoring index for this filter combination
@ -197,7 +197,7 @@ contains
! the end of this select case will handle that one.
if (n < t % scatt_order(j)) then
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) + score
t % results(score_index, filter_index) % value + score
end if
end do
j = j + t % scatt_order(j)
@ -229,7 +229,7 @@ contains
! mu*Sigma_s)).
score = last_wgt / (3.0_8 * score * (material_xs % total - &
mu * score))
mu * score))
case (SCORE_N_1N)
! Skip any event where the particle didn't scatter
@ -299,7 +299,7 @@ contains
! bin. However, in the case of fission, since multiple
! fission neutrons were emitted with different energies,
! multiple outgoing energy bins may have been scored to. The
! following logic treats this special case and scores to
! following logic treats this special case and results to
! multiple bins
call score_fission_eout(t, score_index)
@ -322,13 +322,12 @@ contains
case default
message = "Invalid score type on tally " // to_str(t % id) // "."
write(*,*) score_index, score_bin, l,j, t % n_score_bins
call fatal_error()
end select
! Add score to tally
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) + score
t % results(score_index, filter_index) % value + score
end do SCORE_LOOP
@ -391,13 +390,14 @@ contains
! change outgoing energy bin
t % matching_bins(i) = binary_search(t % filters(i) % real_bins, &
size(t % filters(i) % real_bins), E_out)
size(t % filters(i) % real_bins), E_out)
! determine scoring index
i_filter = sum((t % matching_bins - 1) * t % stride) + 1
! Add score to tally
call add_to_score(t % scores(i_score, i_filter), score)
t % results(i_score, i_filter) % value = &
t % results(i_score, i_filter) % value + score
end do
! reset outgoing energy bin and score index
@ -461,10 +461,10 @@ contains
end if
! =======================================================================
! CALCULATE SCORES AND ACCUMULATE TALLY
! CALCULATE RESULTS AND ACCUMULATE TALLY
! If we have made it here, we have a scoring combination of bins for this
! tally -- now we need to determine where in the scores array we should
! tally -- now we need to determine where in the results array we should
! be accumulating the tally values
! Determine scoring index for this filter combination
@ -472,7 +472,7 @@ contains
if (t % all_nuclides) then
call score_all_nuclides(tracklength_tallies(curr_ptr % data), flux, &
filter_index)
filter_index)
else
NUCLIDE_BIN_LOOP: do k = 1, t % n_nuclide_bins
@ -509,20 +509,20 @@ contains
select case(score_bin)
case (SCORE_TOTAL)
score = micro_xs(i_nuclide) % total * &
atom_density * flux
atom_density * flux
case (SCORE_SCATTER)
score = (micro_xs(i_nuclide) % total - &
micro_xs(i_nuclide) % absorption) * &
atom_density * flux
micro_xs(i_nuclide) % absorption) * &
atom_density * flux
case (SCORE_ABSORPTION)
score = micro_xs(i_nuclide) % absorption * &
atom_density * flux
atom_density * flux
case (SCORE_FISSION)
score = micro_xs(i_nuclide) % fission * &
atom_density * flux
atom_density * flux
case (SCORE_NU_FISSION)
score = micro_xs(i_nuclide) % nu_fission * &
atom_density * flux
atom_density * flux
case (SCORE_EVENTS)
score = ONE
case default
@ -557,7 +557,8 @@ contains
score_index = (k - 1)*t % n_score_bins + j
! Add score to tally
call add_to_score(t % scores(score_index, filter_index), score)
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
end do SCORE_LOOP
@ -632,7 +633,7 @@ contains
score = micro_xs(i_nuclide) % total * atom_density * flux
case (SCORE_SCATTER)
score = (micro_xs(i_nuclide) % total - &
micro_xs(i_nuclide) % absorption) * atom_density * flux
micro_xs(i_nuclide) % absorption) * atom_density * flux
case (SCORE_ABSORPTION)
score = micro_xs(i_nuclide) % absorption * atom_density * flux
case (SCORE_FISSION)
@ -651,7 +652,8 @@ contains
score_index = (i_nuclide - 1)*t % n_score_bins + j
! Add score to tally
call add_to_score(t % scores(score_index, filter_index), score)
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
end do SCORE_LOOP
@ -691,7 +693,8 @@ contains
score_index = n_nuclides_total*t % n_score_bins + j
! Add score to tally
call add_to_score(t % scores(score_index, filter_index), score)
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
end do MATERIAL_SCORE_LOOP
@ -775,27 +778,27 @@ contains
! determine next universe bin
! TODO: Account for multiple universes when performing this filter
t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
p % coord % universe, i_tally)
p % coord % universe, i_tally)
case (FILTER_MATERIAL)
t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, &
p % material, i_tally)
p % material, i_tally)
case (FILTER_CELL)
! determine next cell bin
! TODO: Account for cells in multiple levels when performing this filter
t % matching_bins(i) = get_next_bin(FILTER_CELL, &
p % coord % cell, i_tally)
p % coord % cell, i_tally)
case (FILTER_CELLBORN)
! determine next cellborn bin
t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, &
p % cell_born, i_tally)
p % cell_born, i_tally)
case (FILTER_SURFACE)
! determine next surface bin
t % matching_bins(i) = get_next_bin(FILTER_SURFACE, &
p % surface, i_tally)
p % surface, i_tally)
case (FILTER_ENERGYIN)
! determine incoming energy bin
@ -803,12 +806,12 @@ contains
! check if energy of the particle is within energy bins
if (p % E < t % filters(i) % real_bins(1) .or. &
p % E > t % filters(i) % real_bins(k + 1)) then
p % E > t % filters(i) % real_bins(k + 1)) then
t % matching_bins(i) = NO_BIN_FOUND
else
! search to find incoming energy bin
t % matching_bins(i) = binary_search(t % filters(i) % real_bins, &
k + 1, p % E)
k + 1, p % E)
end if
end select
@ -929,25 +932,25 @@ contains
select case(score_bin)
case (SCORE_TOTAL)
score = micro_xs(i_nuclide) % total * &
atom_density * flux
atom_density * flux
case (SCORE_SCATTER)
score = (micro_xs(i_nuclide) % total - &
micro_xs(i_nuclide) % absorption) * &
atom_density * flux
micro_xs(i_nuclide) % absorption) * &
atom_density * flux
case (SCORE_ABSORPTION)
score = micro_xs(i_nuclide) % absorption * &
atom_density * flux
atom_density * flux
case (SCORE_FISSION)
score = micro_xs(i_nuclide) % fission * &
atom_density * flux
atom_density * flux
case (SCORE_NU_FISSION)
score = micro_xs(i_nuclide) % nu_fission * &
atom_density * flux
atom_density * flux
case (SCORE_EVENTS)
score = ONE
case default
message = "Invalid score type on tally " // &
to_str(t % id) // "."
to_str(t % id) // "."
call fatal_error()
end select
@ -970,7 +973,7 @@ contains
score = ONE
case default
message = "Invalid score type on tally " // &
to_str(t % id) // "."
to_str(t % id) // "."
call fatal_error()
end select
end if
@ -979,7 +982,8 @@ contains
score_index = (b - 1)*t % n_score_bins + j
! Add score to tally
call add_to_score(t % scores(score_index, filter_index), score)
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
end do SCORE_LOOP
@ -1028,27 +1032,27 @@ contains
! determine next universe bin
! TODO: Account for multiple universes when performing this filter
t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
p % coord % universe, i_tally)
p % coord % universe, i_tally)
case (FILTER_MATERIAL)
t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, &
p % material, i_tally)
p % material, i_tally)
case (FILTER_CELL)
! determine next cell bin
! TODO: Account for cells in multiple levels when performing this filter
t % matching_bins(i) = get_next_bin(FILTER_CELL, &
p % coord % cell, i_tally)
p % coord % cell, i_tally)
case (FILTER_CELLBORN)
! determine next cellborn bin
t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, &
p % cell_born, i_tally)
p % cell_born, i_tally)
case (FILTER_SURFACE)
! determine next surface bin
t % matching_bins(i) = get_next_bin(FILTER_SURFACE, &
p % surface, i_tally)
p % surface, i_tally)
case (FILTER_ENERGYIN)
! determine incoming energy bin
@ -1063,12 +1067,12 @@ contains
! check if energy of the particle is within energy bins
if (E < t % filters(i) % real_bins(1) .or. &
E > t % filters(i) % real_bins(n + 1)) then
E > t % filters(i) % real_bins(n + 1)) then
t % matching_bins(i) = NO_BIN_FOUND
else
! search to find incoming energy bin
t % matching_bins(i) = binary_search(t % filters(i) % real_bins, &
n + 1, E)
n + 1, E)
end if
case (FILTER_ENERGYOUT)
@ -1077,12 +1081,12 @@ contains
! check if energy of the particle is within energy bins
if (E < t % filters(i) % real_bins(1) .or. &
E > t % filters(i) % real_bins(n + 1)) then
E > t % filters(i) % real_bins(n + 1)) then
t % matching_bins(i) = NO_BIN_FOUND
else
! search to find incoming energy bin
t % matching_bins(i) = binary_search(t % filters(i) % real_bins, &
n + 1, p % E)
n + 1, p % E)
end if
end select
@ -1173,14 +1177,14 @@ contains
n = t % filters(j) % n_bins
! check if energy of the particle is within energy bins
if (p % E < t % filters(j) % real_bins(1) .or. &
p % E > t % filters(j) % real_bins(n + 1)) then
p % E > t % filters(j) % real_bins(n + 1)) then
curr_ptr => curr_ptr % next ! select next tally
cycle
end if
! search to find incoming energy bin
t % matching_bins(j) = binary_search(t % filters(j) % real_bins, &
n + 1, p % E)
n + 1, p % E)
end if
! =======================================================================
@ -1198,9 +1202,10 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_TOP
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
call add_to_score(t % scores(1, filter_index), p % wgt)
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
end if
end do
else
@ -1209,9 +1214,10 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_TOP
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
call add_to_score(t % scores(1, filter_index), p % wgt)
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
end if
end do
end if
@ -1225,9 +1231,10 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_FRONT
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
call add_to_score(t % scores(1, filter_index), p % wgt)
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
end if
end do
else
@ -1236,9 +1243,10 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_FRONT
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
call add_to_score(t % scores(1, filter_index), p % wgt)
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
end if
end do
end if
@ -1252,9 +1260,10 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_RIGHT
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
call add_to_score(t % scores(1, filter_index), p % wgt)
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
end if
end do
else
@ -1263,9 +1272,10 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_RIGHT
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
call add_to_score(t % scores(1, filter_index), p % wgt)
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
end if
end do
end if
@ -1316,7 +1326,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_RIGHT
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
ijk0(1) = ijk0(1) + 1
xyz_cross(1) = xyz_cross(1) + m % width(1)
@ -1328,7 +1338,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_RIGHT
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
end if
elseif (distance == d(2)) then
@ -1338,7 +1348,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_FRONT
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
ijk0(2) = ijk0(2) + 1
xyz_cross(2) = xyz_cross(2) + m % width(2)
@ -1350,7 +1360,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_FRONT
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
end if
else if (distance == d(3)) then
@ -1360,7 +1370,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_TOP
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
ijk0(3) = ijk0(3) + 1
xyz_cross(3) = xyz_cross(3) + m % width(3)
@ -1372,7 +1382,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_TOP
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
end if
end if
@ -1383,13 +1393,14 @@ contains
! Check for errors
if (filter_index <= 0 .or. filter_index > &
t % total_filter_bins) then
t % total_filter_bins) then
message = "Score index outside range."
call fatal_error()
end if
! Add to surface current tally
call add_to_score(t % scores(1, filter_index), p % wgt)
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
end if
! Calculate new coordinates
@ -1440,7 +1451,7 @@ contains
end if
i_tally_check = tally_maps(filter_type) % items(filter_value) % &
elements(position(filter_type)) % index_tally
elements(position(filter_type)) % index_tally
if (i_tally_check > i_tally) then
! Since the index being checked against is greater than the index we
@ -1452,7 +1463,7 @@ contains
elseif (i_tally_check == i_tally) then
! Found a match
bin = tally_maps(filter_type) % items(filter_value) % &
elements(position(filter_type)) % index_bin
elements(position(filter_type)) % index_bin
return
end if
@ -1484,21 +1495,23 @@ contains
! Accumulate on master only unless run is not reduced then do it on all
if (master .or. (.not. reduce_tallies)) then
! Accumulate scores for each tally
! Accumulate results for each tally
curr_ptr => active_tallies
do while(associated(curr_ptr))
call accumulate_tally(tallies(curr_ptr % data))
curr_ptr => curr_ptr % next
end do
if (run_mode == MODE_CRITICALITY) then
! Before accumulating scores for global_tallies, we need to get the
! current batch estimate of k_analog for displaying to output
if (active_batches) k_batch(current_batch) = global_tallies(K_ANALOG) % value
if (run_mode == MODE_EIGENVALUE) then
! Get the current batch estimate of k_analog for displaying to output
! --- this has to be performed after reduce_tally_values and before
! accumulate_result
k_batch(current_batch) = global_tallies(K_ANALOG) % value
end if
! Accumulate scores for global tallies
if (active_batches) call accumulate_score(global_tallies)
! Accumulate results for global tallies
if (active_batches) call accumulate_result(global_tallies)
end if
if (associated(curr_ptr)) nullify(curr_ptr)
@ -1515,7 +1528,7 @@ contains
integer :: n ! number of filter bins
integer :: m ! number of score bins
integer :: n_bins ! total number of bins
real(8), allocatable :: tally_temp(:,:) ! contiguous array of scores
real(8), allocatable :: tally_temp(:,:) ! contiguous array of results
real(8) :: global_temp(N_GLOBAL_TALLIES)
real(8) :: dummy ! temporary receive buffer for non-root reduces
type(TallyObject), pointer :: t => null()
@ -1531,23 +1544,23 @@ contains
allocate(tally_temp(m,n))
tally_temp = t % scores(:,:) % value
tally_temp = t % results(:,:) % value
if (master) then
! 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)
MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
! Transfer values to value on master
t % scores(:,:) % value = tally_temp
t % results(:,:) % value = tally_temp
else
! Receive buffer not significant at other processors
call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, &
MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
! Reset value on other processors
t % scores(:,:) % value = 0
t % results(:,:) % value = 0
end if
deallocate(tally_temp)
@ -1560,14 +1573,14 @@ contains
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, global_temp, N_GLOBAL_TALLIES, &
MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
! Transfer values back to global_tallies on master
global_tallies(:) % value = global_temp
else
! Receive buffer not significant at other processors
call MPI_REDUCE(global_temp, dummy, N_GLOBAL_TALLIES, &
MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
! Reset value on other processors
global_tallies(:) % value = ZERO
@ -1578,11 +1591,11 @@ contains
! last realization
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, total_weight, 1, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
0, MPI_COMM_WORLD, mpi_err)
else
! Receive buffer not significant at other processors
call MPI_REDUCE(total_weight, dummy, 1, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
0, MPI_COMM_WORLD, mpi_err)
end if
if (associated(curr_ptr)) nullify(curr_ptr)
@ -1605,14 +1618,14 @@ contains
t % n_realizations = t % n_realizations + n_procs
end if
! Accumulate each TallyScore
call accumulate_score(t % scores)
! Accumulate each TallyResult
call accumulate_result(t % results)
end subroutine accumulate_tally
!===============================================================================
! TALLY_STATISTICS computes the mean and standard deviation of the mean of each
! tally and stores them in the val and val_sq attributes of the TallyScores
! tally and stores them in the val and val_sq attributes of the TallyResults
! respectively
!===============================================================================
@ -1625,87 +1638,72 @@ contains
do i = 1, n_tallies
t => tallies(i)
call statistics_score(t % scores, t % n_realizations)
call statistics_result(t % results, t % n_realizations)
end do
! Calculate statistics for global tallies
call statistics_score(global_tallies, n_realizations)
call statistics_result(global_tallies, n_realizations)
end subroutine tally_statistics
!===============================================================================
! ADD_TO_SCORE accumulates a scoring contribution to a specific tally bin and
! specific response function. Note that we don't need to add the square of the
! contribution since that is done at the batch level, not the history level
!===============================================================================
subroutine add_to_score(score, val)
type(TallyScore), intent(inout) :: score
real(8), intent(in) :: val
score % value = score % value + val
end subroutine add_to_score
!===============================================================================
! ACCUMULATE_SCORE accumulates scores from many histories (or many generations)
! ACCUMULATE_RESULT accumulates results from many histories (or many generations)
! into a single realization of a random variable.
!===============================================================================
elemental subroutine accumulate_score(score)
elemental subroutine accumulate_result(this)
type(TallyScore), intent(inout) :: score
type(TallyResult), intent(inout) :: this
real(8) :: val
! Add the sum and square of the sum of contributions from a tally score to
! Add the sum and square of the sum of contributions from a tally result to
! the variables sum and sum_sq. This will later allow us to calculate a
! variance on the tallies.
val = score % value/total_weight
score % sum = score % sum + val
score % sum_sq = score % sum_sq + val*val
val = this % value/total_weight
this % sum = this % sum + val
this % sum_sq = this % sum_sq + val*val
! Reset the single batch estimate
score % value = ZERO
this % value = ZERO
end subroutine accumulate_score
end subroutine accumulate_result
!===============================================================================
! STATISTICS_SCORE determines the sample mean and the standard deviation of the
! mean for a TallyScore.
! STATISTICS_RESULT determines the sample mean and the standard deviation of the
! mean for a TallyResult.
!===============================================================================
elemental subroutine statistics_score(score, n)
elemental subroutine statistics_result(this, n)
type(TallyScore), intent(inout) :: score
integer, intent(in) :: n
type(TallyResult), intent(inout) :: this
integer, intent(in) :: n
! Calculate sample mean and standard deviation of the mean -- note that we
! have used Bessel's correction so that the estimator of the variance of the
! sample mean is unbiased.
score % sum = score % sum/n
score % sum_sq = sqrt((score % sum_sq/n - score % sum * &
score % sum) / (n - 1))
this % sum = this % sum/n
this % sum_sq = sqrt((this % sum_sq/n - this % sum * &
this % sum) / (n - 1))
end subroutine statistics_score
end subroutine statistics_result
!===============================================================================
! RESET_SCORE zeroes out the value and accumulated sum and sum-squared for a
! single TallyScore.
! RESET_RESULT zeroes out the value and accumulated sum and sum-squared for a
! single TallyResult.
!===============================================================================
elemental subroutine reset_score(score)
elemental subroutine reset_result(this)
type(TallyScore), intent(inout) :: score
type(TallyResult), intent(inout) :: this
score % value = ZERO
score % sum = ZERO
score % sum_sq = ZERO
this % value = ZERO
this % sum = ZERO
this % sum_sq = ZERO
end subroutine reset_score
end subroutine reset_result
!===============================================================================
! SETUP_ACTIVE_USERTALLIES
@ -1909,7 +1907,7 @@ contains
end if
do i = n_cmfd_current_tallies + n_user_current_tallies, &
n_user_current_tallies + 1, -1
n_user_current_tallies + 1, -1
! allocate node
allocate(curr_ptr)

20
src/tally_header.F90
View file

@ -35,14 +35,14 @@ module tally_header
end type TallyMap
!===============================================================================
! TALLYSCORE provides accumulation of scores in a particular tally bin
! TALLYRESULT provides accumulation of results in a particular tally bin
!===============================================================================
type TallyScore
type TallyResult
real(8) :: value = 0.
real(8) :: sum = 0.
real(8) :: sum_sq = 0.
end type TallyScore
end type TallyResult
!===============================================================================
! TALLYFILTER describes a filter that limits what events score to a tally. For
@ -59,8 +59,8 @@ module tally_header
!===============================================================================
! TALLYOBJECT describes a user-specified tally. The region of phase space to
! tally in is given by the TallyBins and the scores are stored in a TallyScore
! array.
! tally in is given by the TallyFilters and the results are stored in a
! TallyResult array.
!===============================================================================
type TallyObject
@ -77,10 +77,10 @@ module tally_header
integer :: n_filters ! Number of filters
type(TallyFilter), allocatable :: filters(:) ! Filter data (type/bins)
! The stride attribute is used for determining the index in the scores
! The stride attribute is used for determining the index in the results
! array for a matching_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
! mapped onto one dimension in the results array, the stride attribute gives
! the stride for a given filter type within the results array
integer, allocatable :: matching_bins(:)
integer, allocatable :: stride(:)
@ -105,14 +105,14 @@ module tally_header
integer, allocatable :: scatt_order(:)
integer :: n_nonPN_score_bins = 0
! Scores for each bin -- the first dimenion of the array is for scores
! Results for each bin -- the first dimension of the array is for scores
! (e.g. flux, total reaction rate, fission reaction rate, etc.) and the
! second dimension of the array is for the combination of filters
! (e.g. specific cell, specific energy group, etc.)
integer :: total_filter_bins
integer :: total_score_bins
type(TallyScore), allocatable :: scores(:,:)
type(TallyResult), allocatable :: results(:,:)
! reset property - allows a tally to be reset after every batch
logical :: reset = .false.

8
src/tally_initialize.F90
View file

@ -25,7 +25,7 @@ contains
!===============================================================================
! SETUP_TALLY_ARRAYS allocates and populates several member arrays of the
! TallyObject derived type, including stride, matching_bins, and scores.
! TallyObject derived type, including stride, matching_bins, and results.
!===============================================================================
subroutine setup_tally_arrays()
@ -57,8 +57,8 @@ contains
t % total_filter_bins = n
t % total_score_bins = t % n_score_bins * t % n_nuclide_bins
! Allocate scores array
allocate(t % scores(t % total_score_bins, t % total_filter_bins))
! Allocate results array
allocate(t % results(t % total_score_bins, t % total_filter_bins))
end do TALLY_LOOP
@ -68,7 +68,7 @@ contains
! SETUP_TALLY_MAPS 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
! subroutine also sets the stride attribute for each tally as well as allocating
! storage for the scores array.
! storage for the results array.
!===============================================================================
subroutine setup_tally_maps()

2
src/templates/settings.rnc
View file

@ -2,7 +2,7 @@ element settings {
element confidence_intervals { ( "on" | "off" ) }? &
(
element criticality {
element eigenvalue {
(element batches { xsd:positiveInteger } |
attribute batches { xsd:positiveInteger }) &
(element inactive { xsd:nonNegativeInteger } |

5
src/templates/settings_t.xml
View file

@ -42,9 +42,9 @@
</typedef>
<variable name="confidence_intervals_" tag="confidence_intervals" type="word" length="3" />
<variable name="criticality_" tag="criticality" type="run_parameters_xml" />
<variable name="cross_sections_" tag="cross_sections" type="word" length="255" />
<variable name="cutoff_" tag="cutoff" type="cutoff_xml" dimension="1" />
<variable name="eigenvalue_" tag="eigenvalue" type="run_parameters_xml" />
<variable name="energy_grid_" tag="energy_grid" type="word" length="7" />
<variable name="entropy_" tag="entropy" type="mesh_xml" dimension="1" />
<variable name="fixed_source_" tag="fixed_source" type="run_parameters_xml" />
@ -60,4 +60,7 @@
<variable name="uniform_fs_" tag="uniform_fs" type="mesh_xml" dimension="1" />
<variable name="verbosity_" tag="verbosity" type="integer" />
<!-- Check for old criticality tag -->
<variable name="criticality_" tag="criticality" type="run_parameters_xml" />
</template>

43
src/utils/statepoint.py
View file

@ -12,9 +12,10 @@ filter_types = {1: 'universe', 2: 'material', 3: 'cell', 4: 'cellborn',
score_types = {-1: 'flux', -2: 'total', -3: 'scatter', -4: 'nu-scatter',
-5: 'scatter-1', -6: 'scatter-2', -7: 'scatter-3',
-8: 'transport', -9: 'diffusion', -10: 'n1n', -11: 'n2n',
-12: 'n3n', -13: 'n4n', -14: 'absorption', -15: 'fission',
-16: 'nu-fission', -17: 'current', -18: 'events'}
-8: 'scatter-4', -9: 'scatter-5', -10: 'transport',
-11: 'diffusion', -12: 'n1n', -13: 'n2n', -14: 'n3n',
-15: 'n4n', -16: 'absorption', -17: 'fission',
-18: 'nu-fission', -19: 'current', -20: 'events'}
class BinaryFile(object):
def __init__(self, filename):
@ -80,7 +81,7 @@ class StatePoint(BinaryFile):
# Set flags for what data was read
self._metadata = False
self._values = False
self._results = False
self._source = False
# Initialize arrays for meshes and tallies
@ -193,7 +194,7 @@ class StatePoint(BinaryFile):
# Set flag indicating metadata has already been read
self._metadata = True
def read_values(self):
def read_results(self):
# Check whether metadata has been read
if not self._metadata:
self._read_metadata()
@ -214,16 +215,16 @@ class StatePoint(BinaryFile):
for t in self.tallies:
n = t.total_score_bins * t.total_filter_bins
t.values = np.array(self._get_double(2*n))
t.values.shape = (t.total_filter_bins, t.total_score_bins, 2)
t.results = np.array(self._get_double(2*n))
t.results.shape = (t.total_filter_bins, t.total_score_bins, 2)
# Indicate that tally values have been read
self._values = True
# Indicate that tally results have been read
self._results = True
def read_source(self):
# Check whether tally values have been read
if not self._values:
self.read_values()
# Check whether tally results have been read
if not self._results:
self.read_results()
for i in range(self.n_particles):
s = SourceSite()
@ -273,18 +274,18 @@ class StatePoint(BinaryFile):
# Regular tallies
for t in self.tallies:
for i in range(t.values.shape[0]):
for j in range(t.values.shape[1]):
for i in range(t.results.shape[0]):
for j in range(t.results.shape[1]):
# Get sum and sum of squares
s, s2 = t.values[i,j]
s, s2 = t.results[i,j]
# Calculate sample mean and replace value
s /= n
t.values[i,j,0] = s
t.results[i,j,0] = s
# Calculate standard deviation
if s != 0.0:
t.values[i,j,1] = t_value*sqrt((s2/n - s*s)/(n-1))
t.results[i,j,1] = t_value*sqrt((s2/n - s*s)/(n-1))
def get_value(self, tally_index, spec_list, score_index):
"""Returns a tally score given a list of filters to satisfy.
@ -304,14 +305,14 @@ class StatePoint(BinaryFile):
score_index : int
Index corresponding to score for tally, i.e. the second index in
Tally.values[:,:,:].
Tally.results[:,:,:].
"""
# Get Tally object given the index
t = self.tallies[tally_index]
# Initialize index for filter in Tally.values[:,:,:]
# Initialize index for filter in Tally.results[:,:,:]
filter_index = 0
# Loop over specified filters in spec_list
@ -334,7 +335,7 @@ class StatePoint(BinaryFile):
else:
filter_index += f_index*t.filters[f_type].stride
# Return the desired result from Tally.values. This could be the sum and
# Return the desired result from Tally.results. This could be the sum and
# sum of squares, or it could be mean and stdev if self.generate_stdev()
# has been called already.
return t.values[filter_index, score_index]
return t.results[filter_index, score_index]

6
src/utils/statepoint_cmp.py
View file

@ -17,8 +17,8 @@ sp1 = StatePoint(path1)
sp2 = StatePoint(path2)
# Read tally results
sp1.read_values()
sp2.read_values()
sp1.read_results()
sp2.read_results()
# Compare header information
assert sp1.revision == sp2.revision
@ -71,7 +71,7 @@ for t1, t2 in zip(sp1.tallies, sp2.tallies):
assert t1.scores == t2.scores
# Compare tally results
assert_allclose(t1.values, t2.values)
assert_allclose(t1.results, t2.results)
# If criticality, compare source sites
if sp1.run_mode == 2: