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

Use multidimensional array for tally results rather than derived type

Browse source
This commit is contained in:
Paul Romano 2016-11-14 11:26:08 -06:00
parent 288521a217
commit fb009c6b5e
24 changed files with 331 additions and 502 deletions
Download patch file Download diff file Expand all files Collapse all files

10
docs/source/io_formats/statepoint.rst
View file

@ -224,12 +224,12 @@ if run_mode == 'k-eigenvalue':
Tallying moment orders for Legendre and spherical harmonic tally expansions
(*e.g.*, 'P2', 'Y1,2', etc.).
**/tallies/tally <uid>/results** (Compound type)
**/tallies/tally <uid>/results** (*double[][][2]*)
Accumulated sum and sum-of-squares for each bin of the i-th tally. This is a
two-dimensional array, the first dimension of which represents combinations
of filter bins and the second dimensions of which represents scoring
bins. Each element of the array has fields 'sum' and 'sum_sq'.
Accumulated sum and sum-of-squares for each bin of the i-th tally. The first
dimension represents combinations of filter bins, the second dimensions
represents scoring bins, and the third dimension has two entries for the sum
and the sum-of-squares.
**/source_present** (*int*)

8
openmc/particle_restart.py
View file

@ -1,3 +1,5 @@
import h5py
class Particle(object):
"""Information used to restart a specific particle that caused a simulation to
fail.
@ -33,12 +35,6 @@ class Particle(object):
"""
def __init__(self, filename):
import h5py
if h5py.__version__ == '2.6.0':
raise ImportError("h5py 2.6.0 has a known bug which makes it "
"incompatible with OpenMC's HDF5 files. "
"Please switch to a different version.")
self._f = h5py.File(filename, 'r')
# Ensure filetype and revision are correct

13
openmc/statepoint.py
View file

@ -5,6 +5,7 @@ import warnings
import glob
import numpy as np
import h5py
import openmc
import openmc.checkvalue as cv
@ -104,12 +105,6 @@ class StatePoint(object):
"""
def __init__(self, filename, autolink=True):
import h5py
if h5py.__version__ == '2.6.0':
raise ImportError("h5py 2.6.0 has a known bug which makes it "
"incompatible with OpenMC's HDF5 files. "
"Please switch to a different version.")
self._f = h5py.File(filename, 'r')
# Ensure filetype and revision are correct
@ -209,13 +204,13 @@ class StatePoint(object):
def global_tallies(self):
if self._global_tallies is None:
data = self._f['global_tallies'].value
gt = np.zeros_like(data, dtype=[
gt = np.zeros(data.shape[0], dtype=[
('name', 'a14'), ('sum', 'f8'), ('sum_sq', 'f8'),
('mean', 'f8'), ('std_dev', 'f8')])
gt['name'] = ['k-collision', 'k-absorption', 'k-tracklength',
'leakage']
gt['sum'] = data['sum']
gt['sum_sq'] = data['sum_sq']
gt['sum'] = data[:,1]
gt['sum_sq'] = data[:,2]
# Calculate mean and sample standard deviation of mean
n = self.n_realizations

10
openmc/summary.py
View file

@ -2,6 +2,7 @@ from collections import Iterable
import re
import numpy as np
import h5py
import openmc
from openmc.region import Region
@ -23,15 +24,6 @@ class Summary(object):
"""
def __init__(self, filename):
# A user may not have h5py, but they can still use the rest of the
# Python API so we'll only try to import h5py if the user actually inits
# a Summary object.
import h5py
if h5py.__version__ == '2.6.0':
raise ImportError("h5py 2.6.0 has a known bug which makes it "
"incompatible with OpenMC's HDF5 files. "
"Please switch to a different version.")
openmc.reset_auto_ids()
if not filename.endswith(('.h5', '.hdf5')):

13
openmc/tallies.py
View file

@ -13,6 +13,7 @@ from xml.etree import ElementTree as ET
from six import string_types
import numpy as np
import h5py
import openmc
import openmc.checkvalue as cv
@ -269,20 +270,14 @@ class Tally(object):
return None
if not self._results_read:
import h5py
if h5py.__version__ == '2.6.0':
raise ImportError("h5py 2.6.0 has a known bug which makes it "
"incompatible with OpenMC's HDF5 files. "
"Please switch to a different version.")
# Open the HDF5 statepoint file
f = h5py.File(self._sp_filename, 'r')
# Extract Tally data from the file
data = f['tallies/tally {0}/results'.format(
self.id)].value
sum = data['sum']
sum_sq = data['sum_sq']
sum = data[:,:,0]
sum_sq = data[:,:,1]
# Reshape the results arrays
sum = np.reshape(sum, self.shape)
@ -1726,8 +1721,6 @@ class Tally(object):
# HDF5 binary file
if format == 'hdf5':
import h5py
filename = directory + '/' + filename + '.h5'
if append:

3
openmc/volume.py
View file

@ -5,6 +5,7 @@ from warnings import warn
import numpy as np
import pandas as pd
import h5py
import openmc
import openmc.checkvalue as cv
@ -187,8 +188,6 @@ class VolumeCalculation(object):
Results of the stochastic volume calculation
"""
import h5py
with h5py.File(filename, 'r') as f:
domain_type = f.attrs['domain_type'].decode()
samples = f.attrs['samples']

39
src/cmfd_data.F90
View file

@ -164,7 +164,7 @@ contains
* t%stride) + 1
! Get flux
flux = t % results(1,score_index) % sum
flux = t % results(RESULT_SUM,1,score_index)
cmfd % flux(h,i,j,k) = flux
! Detect zero flux, abort if located
@ -175,10 +175,10 @@ contains
end if
! Get total rr and convert to total xs
cmfd % totalxs(h,i,j,k) = t % results(2,score_index) % sum / flux
cmfd % totalxs(h,i,j,k) = t % results(RESULT_SUM,2,score_index) / flux
! Get p1 scatter rr and convert to p1 scatter xs
cmfd % p1scattxs(h,i,j,k) = t % results(3,score_index) % sum / flux
cmfd % p1scattxs(h,i,j,k) = t % results(RESULT_SUM,3,score_index) / flux
! Calculate diffusion coefficient
cmfd % diffcof(h,i,j,k) = ONE/(3.0_8*(cmfd % totalxs(h,i,j,k) - &
@ -211,19 +211,18 @@ contains
* t%stride) + 1
! Get scattering
cmfd % scattxs(h,g,i,j,k) = t % results(1,score_index) % sum /&
cmfd % scattxs(h,g,i,j,k) = t % results(RESULT_SUM,1,score_index) /&
cmfd % flux(h,i,j,k)
! Get nu-fission
cmfd % nfissxs(h,g,i,j,k) = t % results(2,score_index) % sum /&
cmfd % nfissxs(h,g,i,j,k) = t % results(RESULT_SUM,2,score_index) /&
cmfd % flux(h,i,j,k)
! Bank source
cmfd % openmc_src(g,i,j,k) = cmfd % openmc_src(g,i,j,k) + &
t % results(2,score_index) % sum
t % results(RESULT_SUM,2,score_index)
cmfd % keff_bal = cmfd % keff_bal + &
t % results(2,score_index) % sum / &
dble(t % n_realizations)
t % results(RESULT_SUM,2,score_index) / t % n_realizations
end do INGROUP
@ -243,67 +242,67 @@ contains
matching_bins(i_filter_surf) = OUT_LEFT
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
cmfd % current(1,h,i,j,k) = t % results(1,score_index) % sum
cmfd % current(1,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
matching_bins(i_filter_surf) = IN_LEFT
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
cmfd % current(2,h,i,j,k) = t % results(1,score_index) % sum
cmfd % current(2,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
! Right surface
matching_bins(i_filter_surf) = IN_RIGHT
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
cmfd % current(3,h,i,j,k) = t % results(1,score_index) % sum
cmfd % current(3,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
matching_bins(i_filter_surf) = OUT_RIGHT
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
cmfd % current(4,h,i,j,k) = t % results(1,score_index) % sum
cmfd % current(4,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
! Back surface
matching_bins(i_filter_surf) = OUT_BACK
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
cmfd % current(5,h,i,j,k) = t % results(1,score_index) % sum
cmfd % current(5,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
matching_bins(i_filter_surf) = IN_BACK
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
cmfd % current(6,h,i,j,k) = t % results(1,score_index) % sum
cmfd % current(6,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
! Front surface
matching_bins(i_filter_surf) = IN_FRONT
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
cmfd % current(7,h,i,j,k) = t % results(1,score_index) % sum
cmfd % current(7,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
matching_bins(i_filter_surf) = OUT_FRONT
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
cmfd % current(8,h,i,j,k) = t % results(1,score_index) % sum
cmfd % current(8,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
! Bottom surface
matching_bins(i_filter_surf) = OUT_BOTTOM
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
cmfd % current(9,h,i,j,k) = t % results(1,score_index) % sum
cmfd % current(9,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
matching_bins(i_filter_surf) = IN_BOTTOM
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
cmfd % current(10,h,i,j,k) = t % results(1,score_index) % sum
cmfd % current(10,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
! Top surface
matching_bins(i_filter_surf) = IN_TOP
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
cmfd % current(11,h,i,j,k) = t % results(1,score_index) % sum
cmfd % current(11,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
matching_bins(i_filter_surf) = OUT_TOP
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
cmfd % current(12,h,i,j,k) = t % results(1,score_index) % sum
cmfd % current(12,h,i,j,k) = t % results(RESULT_SUM,1,score_index)
end if TALLY

12
src/cmfd_execute.F90
View file

@ -365,22 +365,18 @@ contains
subroutine cmfd_tally_reset()
use global, only: n_cmfd_tallies, cmfd_tallies
use global, only: cmfd_tallies
use output, only: write_message
use tally, only: reset_result
integer :: i ! loop counter
! Print message
call write_message("CMFD tallies reset", 7)
! Begin loop around CMFD tallies
do i = 1, n_cmfd_tallies
! Reset that tally
! Reset CMFD tallies
do i = 1, size(cmfd_tallies)
cmfd_tallies(i) % n_realizations = 0
call reset_result(cmfd_tallies(i) % results)
cmfd_tallies(i) % results(:,:,:) = ZERO
end do
end subroutine cmfd_tally_reset

6
src/constants.F90
View file

@ -269,6 +269,12 @@ module constants
! ============================================================================
! TALLY-RELATED CONSTANTS
! Tally result entries
integer, parameter :: &
RESULT_VALUE = 1, &
RESULT_SUM = 2, &
RESULT_SUM_SQ = 3
! Tally type
integer, parameter :: &
TALLY_VOLUME = 1, &

14
src/eigenvalue.F90
View file

@ -374,7 +374,7 @@ contains
subroutine calculate_generation_keff()
! Get keff for this generation by subtracting off the starting value
keff_generation = global_tallies(K_TRACKLENGTH) % value - keff_generation
keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH) - keff_generation
#ifdef MPI
! Combine values across all processors
@ -466,14 +466,14 @@ contains
k_combined = ZERO
! Copy estimates of k-effective and its variance (not variance of the mean)
kv(1) = global_tallies(K_COLLISION) % sum / n
kv(2) = global_tallies(K_ABSORPTION) % sum / n
kv(3) = global_tallies(K_TRACKLENGTH) % sum / n
cov(1,1) = (global_tallies(K_COLLISION) % sum_sq - &
kv(1) = global_tallies(RESULT_SUM, K_COLLISION) / n
kv(2) = global_tallies(RESULT_SUM, K_ABSORPTION) / n
kv(3) = global_tallies(RESULT_SUM, K_TRACKLENGTH) / n
cov(1,1) = (global_tallies(RESULT_SUM_SQ, K_COLLISION) - &
n * kv(1) * kv(1)) / (n - 1)
cov(2,2) = (global_tallies(K_ABSORPTION) % sum_sq - &
cov(2,2) = (global_tallies(RESULT_SUM_SQ, K_ABSORPTION) - &
n * kv(2) * kv(2)) / (n - 1)
cov(3,3) = (global_tallies(K_TRACKLENGTH) % sum_sq - &
cov(3,3) = (global_tallies(RESULT_SUM_SQ, K_TRACKLENGTH) - &
n * kv(3) * kv(3)) / (n - 1)
! Calculate covariances based on sums with Bessel's correction

4
src/finalize.F90
View file

@ -9,7 +9,7 @@ module finalize
use message_passing
#endif
use hdf5_interface, only: hdf5_bank_t, hdf5_tallyresult_t
use hdf5_interface, only: hdf5_bank_t
use hdf5, only: h5tclose_f, h5close_f
implicit none
@ -53,7 +53,6 @@ contains
call free_memory()
! Release compound datatypes
call h5tclose_f(hdf5_tallyresult_t, hdf5_err)
call h5tclose_f(hdf5_bank_t, hdf5_err)
! Close FORTRAN interface.
@ -62,7 +61,6 @@ contains
#ifdef MPI
! Free all MPI types
call MPI_TYPE_FREE(MPI_BANK, mpi_err)
call MPI_TYPE_FREE(MPI_TALLYRESULT, mpi_err)
! If MPI is in use and enabled, terminate it
call MPI_FINALIZE(mpi_err)

16
src/global.F90
View file

@ -1,5 +1,11 @@
module global
use, intrinsic :: ISO_C_BINDING
#ifdef MPIF08
use mpi_f08
#endif
use bank_header, only: Bank
use cmfd_header
use constants
@ -14,15 +20,11 @@ module global
use set_header, only: SetInt
use surface_header, only: SurfaceContainer
use source_header, only: SourceDistribution
use tally_header, only: TallyObject, TallyResult
use tally_header, only: TallyObject
use trigger_header, only: KTrigger
use timer_header, only: Timer
use volume_header, only: VolumeCalculation
#ifdef MPIF08
use mpi_f08
#endif
implicit none
! ============================================================================
@ -164,7 +166,7 @@ module global
! 3) track-length estimate of k-eff
! 4) leakage fraction
type(TallyResult), allocatable, target :: global_tallies(:)
real(C_DOUBLE), allocatable, target :: global_tallies(:,:)
! It is possible to protect accumulate operations on global tallies by using
! an atomic update. However, when multiple threads accumulate to the same
@ -272,10 +274,8 @@ module global
integer :: mpi_err ! MPI error code
#ifdef MPIF08
type(MPI_Datatype) :: MPI_BANK
type(MPI_Datatype) :: MPI_TALLYRESULT
#else
integer :: MPI_BANK ! MPI datatype for fission bank
integer :: MPI_TALLYRESULT ! MPI datatype for TallyResult
#endif
#ifdef _OPENMP

130
src/hdf5_interface.F90
View file

@ -16,7 +16,6 @@ module hdf5_interface
use h5lt
use error, only: fatal_error
use tally_header, only: TallyResult
#ifdef PHDF5
use message_passing, only: MPI_COMM_WORLD, MPI_INFO_NULL
#endif
@ -24,7 +23,6 @@ module hdf5_interface
implicit none
private
integer(HID_T), public :: hdf5_tallyresult_t ! Compound type for TallyResult
integer(HID_T), public :: hdf5_bank_t ! Compound type for Bank
integer(HID_T), public :: hdf5_integer8_t ! type for integer(8)
@ -42,8 +40,6 @@ module hdf5_interface
module procedure write_long
module procedure write_string
module procedure write_string_1D
module procedure write_tally_result_1D
module procedure write_tally_result_2D
end interface write_dataset
interface read_dataset
@ -60,8 +56,6 @@ module hdf5_interface
module procedure read_long
module procedure read_string
module procedure read_string_1D
module procedure read_tally_result_1D
module procedure read_tally_result_2D
module procedure read_complex_2D
end interface read_dataset
@ -2063,130 +2057,6 @@ contains
call h5ltset_attribute_string_f(group_id, var, attr_type, attr_str, hdf5_err)
end subroutine write_attribute_string
!===============================================================================
! WRITE_TALLY_RESULT writes an OpenMC TallyResult type
!===============================================================================
subroutine write_tally_result_1D(group_id, name, buffer)
integer(HID_T), intent(in) :: group_id
character(*), intent(in) :: name ! name of data
type(TallyResult), intent(in), target :: buffer(:) ! data to write
integer(HSIZE_T) :: dims(1)
dims(:) = shape(buffer)
call write_tally_result_1D_explicit(group_id, dims, name, buffer)
end subroutine write_tally_result_1D
subroutine write_tally_result_1D_explicit(group_id, dims, name, buffer)
integer(HID_T), intent(in) :: group_id
integer(HSIZE_T), intent(in) :: dims(1)
character(*), intent(in) :: name ! name of data
type(TallyResult), intent(in), target :: buffer(dims(1))
integer :: hdf5_err
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
call h5screate_simple_f(1, dims, dspace, hdf5_err)
call h5dcreate_f(group_id, trim(name), hdf5_tallyresult_t, &
dspace, dset, hdf5_err)
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err)
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
end subroutine write_tally_result_1D_explicit
subroutine write_tally_result_2D(group_id, name, buffer)
integer(HID_T), intent(in) :: group_id
character(*), intent(in) :: name ! name of data
type(TallyResult), intent(in), target :: buffer(:,:) ! data to write
integer(HSIZE_T) :: dims(2)
dims(:) = shape(buffer)
call write_tally_result_2D_explicit(group_id, dims, name, buffer)
end subroutine write_tally_result_2D
subroutine write_tally_result_2D_explicit(group_id, dims, name, buffer)
integer(HID_T), intent(in) :: group_id
integer(HSIZE_T), intent(in) :: dims(2)
character(*), intent(in) :: name ! name of data
type(TallyResult), intent(in), target :: buffer(dims(1),dims(2))
integer :: hdf5_err
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
type(c_ptr) :: f_ptr
call h5screate_simple_f(2, dims, dspace, hdf5_err)
call h5dcreate_f(group_id, trim(name), hdf5_tallyresult_t, &
dspace, dset, hdf5_err)
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err)
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
end subroutine write_tally_result_2D_explicit
!===============================================================================
! READ_TALLY_RESULT reads OpenMC TallyResult data
!===============================================================================
subroutine read_tally_result_1D(group_id, name, buffer)
integer(HID_T), intent(in) :: group_id
character(*), intent(in) :: name ! name of data
type(TallyResult), intent(inout), target :: buffer(:) ! read data here
integer(HSIZE_T) :: dims(1)
dims(:) = shape(buffer)
call read_tally_result_1D_explicit(group_id, dims, name, buffer)
end subroutine read_tally_result_1D
subroutine read_tally_result_1D_explicit(group_id, dims, name, buffer)
integer(HID_T), intent(in) :: group_id
integer(HSIZE_T), intent(in) :: dims(1)
character(*), intent(in) :: name ! name of data
type(TallyResult), intent(inout), target :: buffer(dims(1))
integer :: hdf5_err
integer(HID_T) :: dset ! data set handle
type(c_ptr) :: f_ptr
call h5dopen_f(group_id, trim(name), dset, hdf5_err)
f_ptr = c_loc(buffer)
call h5dread_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err)
call h5dclose_f(dset, hdf5_err)
end subroutine read_tally_result_1D_explicit
subroutine read_tally_result_2D(group_id, name, buffer)
integer(HID_T), intent(in) :: group_id
character(*), intent(in) :: name ! name of data
type(TallyResult), intent(inout), target :: buffer(:,:)
integer(HSIZE_T) :: dims(2)
dims(:) = shape(buffer)
call read_tally_result_2D_explicit(group_id, dims, name, buffer)
end subroutine read_tally_result_2D
subroutine read_tally_result_2D_explicit(group_id, dims, name, buffer)
integer(HID_T), intent(in) :: group_id
integer(HSIZE_T), intent(in) :: dims(2)
character(*), intent(in) :: name ! name of data
type(TallyResult), intent(inout), target :: buffer(dims(1),dims(2))
integer :: hdf5_err
integer(HID_T) :: dset ! data set handle
type(c_ptr) :: f_ptr
call h5dopen_f(group_id, trim(name), dset, hdf5_err)
f_ptr = c_loc(buffer)
call h5dread_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err)
call h5dclose_f(dset, hdf5_err)
end subroutine read_tally_result_2D_explicit
subroutine read_attribute_double(buffer, obj_id, name)
real(8), intent(inout), target :: buffer
integer(HID_T), intent(in) :: obj_id

53
src/initialize.F90
View file

@ -12,7 +12,7 @@ module initialize
&BASE_UNIVERSE
use global
use hdf5_interface, only: file_open, read_dataset, file_close, hdf5_bank_t,&
hdf5_tallyresult_t, hdf5_integer8_t
hdf5_integer8_t
use input_xml, only: read_input_xml, cells_in_univ_dict, read_plots_xml
use material_header, only: Material
use mgxs_data, only: read_mgxs, create_macro_xs
@ -22,7 +22,7 @@ module initialize
use state_point, only: load_state_point
use string, only: to_str, starts_with, ends_with, str_to_int
use summary, only: write_summary
use tally_header, only: TallyObject, TallyResult
use tally_header, only: TallyObject
use tally_initialize,only: configure_tallies
use tally_filter
use tally, only: init_tally_routines
@ -169,21 +169,11 @@ contains
integer :: bank_blocks(5) ! Count for each datatype
#ifdef MPIF08
type(MPI_Datatype) :: bank_types(5)
type(MPI_Datatype) :: result_types(1)
type(MPI_Datatype) :: temp_type
#else
integer :: bank_types(5) ! Datatypes
integer :: result_types(1) ! Datatypes
integer :: temp_type ! temporary derived type
#endif
integer(MPI_ADDRESS_KIND) :: bank_disp(5) ! Displacements
integer :: result_blocks(1) ! Count for each datatype
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(TallyResult) :: tr
! Indicate that MPI is turned on
mpi_enabled = .true.
@ -221,36 +211,6 @@ contains
call MPI_TYPE_CREATE_STRUCT(5, bank_blocks, bank_disp, &
bank_types, MPI_BANK, mpi_err)
call MPI_TYPE_COMMIT(MPI_BANK, mpi_err)
! ==========================================================================
! CREATE MPI_TALLYRESULT TYPE
! Determine displacements for MPI_BANK type
call MPI_GET_ADDRESS(tr%value, result_base_disp, mpi_err)
call MPI_GET_ADDRESS(tr%sum, result_disp(1), mpi_err)
! Adjust displacements
result_disp = result_disp - result_base_disp
! 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 = result_disp(1) + 16
call MPI_TYPE_CREATE_RESIZED(temp_type, lower_bound, extent, &
MPI_TALLYRESULT, mpi_err)
! Commit derived type for tally scores
call MPI_TYPE_COMMIT(MPI_TALLYRESULT, mpi_err)
! Free temporary MPI type
call MPI_TYPE_FREE(temp_type, mpi_err)
end subroutine initialize_mpi
#endif
!===============================================================================
@ -259,7 +219,6 @@ contains
subroutine hdf5_initialize()
type(TallyResult), target :: tmp(2) ! temporary TallyResult
type(Bank), target :: tmpb(2) ! temporary Bank
integer :: hdf5_err
integer(HID_T) :: coordinates_t ! HDF5 type for 3 reals
@ -268,14 +227,6 @@ contains
! Initialize FORTRAN interface.
call h5open_f(hdf5_err)
! Create the compound datatype for TallyResult
call h5tcreate_f(H5T_COMPOUND_F, 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_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
call h5tarray_create_f(H5T_NATIVE_DOUBLE, 1, dims, coordinates_t, hdf5_err)

1
src/input_xml.F90
View file

@ -2039,7 +2039,6 @@ contains
type(Library), allocatable :: libraries(:)
type(VectorReal), allocatable :: nuc_temps(:) ! List of T to read for each nuclide
type(VectorReal), allocatable :: sab_temps(:) ! List of T to read for each S(a,b)
character(MAX_LINE_LEN) :: temp_str
real(8), allocatable :: material_temps(:)
logical :: file_exists
character(MAX_FILE_LEN) :: env_variable

3
src/math.F90
View file

@ -1,8 +1,9 @@
module math
use, intrinsic :: ISO_C_BINDING
use constants
use random_lcg, only: prn
use ISO_C_BINDING
implicit none

95
src/output.F90
View file

@ -611,7 +611,7 @@ contains
t_value = t_percentile(ONE - alpha/TWO, n_realizations - 1)
! Adjust sum_sq
global_tallies(:) % sum_sq = t_value * global_tallies(:) % sum_sq
global_tallies(RESULT_SUM_SQ,:) = t_value * global_tallies(RESULT_SUM_SQ,:)
! Adjust combined estimator
if (n_realizations > 3) then
@ -623,26 +623,26 @@ contains
! write global tallies
if (n_realizations > 1) then
if (run_mode == MODE_EIGENVALUE) then
write(ou,102) "k-effective (Collision)", global_tallies(K_COLLISION) &
% sum, global_tallies(K_COLLISION) % sum_sq
write(ou,102) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) &
% sum, global_tallies(K_TRACKLENGTH) % sum_sq
write(ou,102) "k-effective (Absorption)", global_tallies(K_ABSORPTION) &
% sum, global_tallies(K_ABSORPTION) % sum_sq
write(ou,102) "k-effective (Collision)", global_tallies(RESULT_SUM, &
K_COLLISION), global_tallies(RESULT_SUM_SQ, K_COLLISION)
write(ou,102) "k-effective (Track-length)", global_tallies(RESULT_SUM, &
K_TRACKLENGTH), global_tallies(RESULT_SUM_SQ, K_TRACKLENGTH)
write(ou,102) "k-effective (Absorption)", global_tallies(RESULT_SUM, &
K_ABSORPTION), global_tallies(RESULT_SUM_SQ, K_ABSORPTION)
if (n_realizations > 3) write(ou,102) "Combined k-effective", k_combined
end if
write(ou,102) "Leakage Fraction", global_tallies(LEAKAGE) % sum, &
global_tallies(LEAKAGE) % sum_sq
write(ou,102) "Leakage Fraction", global_tallies(RESULT_SUM, LEAKAGE), &
global_tallies(RESULT_SUM_SQ, LEAKAGE)
else
if (master) call warning("Could not compute uncertainties -- only one &
&active batch simulated!")
if (run_mode == MODE_EIGENVALUE) then
write(ou,103) "k-effective (Collision)", global_tallies(K_COLLISION) % sum
write(ou,103) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) % sum
write(ou,103) "k-effective (Absorption)", global_tallies(K_ABSORPTION) % sum
write(ou,103) "k-effective (Collision)", global_tallies(RESULT_SUM, K_COLLISION)
write(ou,103) "k-effective (Track-length)", global_tallies(RESULT_SUM, K_TRACKLENGTH)
write(ou,103) "k-effective (Absorption)", global_tallies(RESULT_SUM, K_ABSORPTION)
end if
write(ou,103) "Leakage Fraction", global_tallies(LEAKAGE) % sum
write(ou,103) "Leakage Fraction", global_tallies(RESULT_SUM, LEAKAGE)
end if
write(ou,*)
@ -765,7 +765,7 @@ contains
end if
! Multiply uncertainty by t-value
t % results % sum_sq = t_value * t % results % sum_sq
t % results(RESULT_SUM_SQ,:,:) = t_value * t % results(RESULT_SUM_SQ,:,:)
end if
! Write header block
@ -876,8 +876,8 @@ contains
score_names(abs(t % score_bins(k)))
write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index) % sum_sq))
to_str(t % results(RESULT_SUM,score_index,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,score_index,filter_index)))
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
score_index = score_index - 1
do n_order = 0, t % moment_order(k)
@ -886,9 +886,8 @@ contains
score_names(abs(t % score_bins(k)))
write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index) &
% sum_sq))
to_str(t % results(RESULT_SUM,score_index,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,score_index,filter_index)))
end do
k = k + t % moment_order(k)
case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, &
@ -902,9 +901,9 @@ contains
// score_names(abs(t % score_bins(k)))
write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index)&
% sum_sq))
to_str(t % results(RESULT_SUM,score_index,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,score_index,&
filter_index)))
end do
end do
k = k + (t % moment_order(k) + 1)**2 - 1
@ -916,8 +915,8 @@ contains
end if
write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index) % sum_sq))
to_str(t % results(RESULT_SUM,score_index,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,score_index,filter_index)))
end select
end do
indent = indent - 2
@ -1020,16 +1019,16 @@ contains
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current on Left", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
matching_bins(i_filter_surf) = IN_LEFT
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current on Left", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
! Right Surface
matching_bins(i_filter_surf) = OUT_RIGHT
@ -1037,16 +1036,16 @@ contains
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current on Right", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current on Right", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
if (n_dim >= 2) then
@ -1056,16 +1055,16 @@ contains
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current on Back", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
matching_bins(i_filter_surf) = IN_BACK
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current on Back", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
! Front Surface
matching_bins(i_filter_surf) = OUT_FRONT
@ -1073,16 +1072,16 @@ contains
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Net Current on Front", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Net Current on Front", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
end if
if (n_dim == 3) then
@ -1092,16 +1091,16 @@ contains
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current on Bottom", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
matching_bins(i_filter_surf) = IN_BOTTOM
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current on Bottom", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
! Top Surface
matching_bins(i_filter_surf) = OUT_TOP
@ -1109,16 +1108,16 @@ contains
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current on Top", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current on Top", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
to_str(t % results(RESULT_SUM,1,filter_index)), &
trim(to_str(t % results(RESULT_SUM_SQ,1,filter_index)))
end if
end do
end do

8
src/physics_mg.F90
View file

@ -119,16 +119,16 @@ contains
! Score implicit absorption estimate of keff
!$omp atomic
global_tallies(K_ABSORPTION) % value = &
global_tallies(K_ABSORPTION) % value + p % absorb_wgt * &
global_tallies(RESULT_VALUE, K_ABSORPTION) = &
global_tallies(RESULT_VALUE, K_ABSORPTION) + p % absorb_wgt * &
material_xs % nu_fission / material_xs % absorption
else
! See if disappearance reaction happens
if (material_xs % absorption > prn() * material_xs % total) then
! Score absorption estimate of keff
!$omp atomic
global_tallies(K_ABSORPTION) % value = &
global_tallies(K_ABSORPTION) % value + p % wgt * &
global_tallies(RESULT_VALUE, K_ABSORPTION) = &
global_tallies(RESULT_VALUE, K_ABSORPTION) + p % wgt * &
material_xs % nu_fission / material_xs % absorption
p % alive = .false.

31
src/simulation.F90
View file

@ -20,8 +20,7 @@ module simulation
use source, only: initialize_source, sample_external_source
use state_point, only: write_state_point, write_source_point
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies, &
reset_result
use tally, only: synchronize_tallies, setup_active_usertallies
use trigger, only: check_triggers
use tracking, only: transport
use volume_calc, only: run_volume_calculations
@ -220,7 +219,7 @@ contains
if (ufs) call count_source_for_ufs()
! Store current value of tracklength k
keff_generation = global_tallies(K_TRACKLENGTH) % value
keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH)
end if
end subroutine initialize_generation
@ -237,24 +236,24 @@ contains
!$omp parallel
!$omp critical
if (run_mode == MODE_EIGENVALUE) then
global_tallies(K_COLLISION) % value = &
global_tallies(K_COLLISION) % value + global_tally_collision
global_tallies(K_ABSORPTION) % value = &
global_tallies(K_ABSORPTION) % value + global_tally_absorption
global_tallies(K_TRACKLENGTH) % value = &
global_tallies(K_TRACKLENGTH) % value + global_tally_tracklength
global_tallies(RESULT_VALUE, K_COLLISION) = &
global_tallies(RESULT_VALUE, K_COLLISION) + global_tally_collision
global_tallies(RESULT_VALUE, K_ABSORPTION) = &
global_tallies(RESULT_VALUE, K_ABSORPTION) + global_tally_absorption
global_tallies(RESULT_VALUE, K_TRACKLENGTH) = &
global_tallies(RESULT_VALUE, K_TRACKLENGTH) + global_tally_tracklength
end if
global_tallies(LEAKAGE) % value = &
global_tallies(LEAKAGE) % value + global_tally_leakage
global_tallies(RESULT_VALUE, LEAKAGE) = &
global_tallies(RESULT_VALUE, LEAKAGE) + global_tally_leakage
!$omp end critical
! reset private tallies
if (run_mode == MODE_EIGENVALUE) then
global_tally_collision = 0
global_tally_absorption = 0
global_tally_tracklength = 0
global_tally_collision = ZERO
global_tally_absorption = ZERO
global_tally_tracklength = ZERO
end if
global_tally_leakage = 0
global_tally_leakage = ZERO
!$omp end parallel
if (run_mode == MODE_EIGENVALUE) then
@ -302,7 +301,7 @@ contains
! Reset global tally results
if (.not. active_batches) then
call reset_result(global_tallies)
global_tallies(:,:) = ZERO
n_realizations = 0
end if

50
src/state_point.F90
View file

@ -353,7 +353,7 @@ contains
! Write sum and sum_sq for each bin
tally_group = open_group(tallies_group, "tally " &
// to_str(tally % id))
call write_dataset(tally_group, "results", tally % results)
call tally % write_results_hdf5(tally_group)
call close_group(tally_group)
end do TALLY_RESULTS
@ -481,7 +481,7 @@ contains
integer :: n_bins ! total number of bins
integer(HID_T) :: tallies_group, tally_group
real(8), allocatable :: tally_temp(:,:,:) ! contiguous array of results
real(8), target :: global_temp(2,N_GLOBAL_TALLIES)
real(8), target :: global_temp(3,N_GLOBAL_TALLIES)
#ifdef MPI
real(8) :: dummy ! temporary receive buffer for non-root reduces
#endif
@ -489,7 +489,7 @@ contains
type(ElemKeyValueII), pointer :: current
type(ElemKeyValueII), pointer :: next
type(TallyObject), pointer :: tally
type(TallyResult), allocatable :: tallyresult_temp(:,:)
type(TallyObject) :: dummy_tally
! ==========================================================================
! COLLECT AND WRITE GLOBAL TALLIES
@ -505,9 +505,8 @@ contains
end if
! Copy global tallies into temporary array for reducing
n_bins = 2 * N_GLOBAL_TALLIES
global_temp(1,:) = global_tallies(:)%sum
global_temp(2,:) = global_tallies(:)%sum_sq
n_bins = 3 * N_GLOBAL_TALLIES
global_temp(:,:) = global_tallies(:,:)
if (master) then
! The MPI_IN_PLACE specifier allows the master to copy values into a
@ -519,20 +518,11 @@ contains
! Transfer values to value on master
if (current_batch == n_max_batches .or. satisfy_triggers) then
global_tallies(:)%sum = global_temp(1,:)
global_tallies(:)%sum_sq = global_temp(2,:)
global_tallies(:,:) = global_temp(:,:)
end if
! Put reduced value in temporary tally result
allocate(tallyresult_temp(N_GLOBAL_TALLIES, 1))
tallyresult_temp(:,1)%sum = global_temp(1,:)
tallyresult_temp(:,1)%sum_sq = global_temp(2,:)
! Write out global tallies sum and sum_sq
call write_dataset(file_id, "global_tallies", tallyresult_temp)
! Deallocate temporary tally result
deallocate(tallyresult_temp)
call write_dataset(file_id, "global_tallies", global_temp)
else
! Receive buffer not significant at other processors
#ifdef MPI
@ -568,15 +558,15 @@ contains
tally => tallies(i)
! Determine size of tally results array
m = size(tally%results, 1)
n = size(tally%results, 2)
m = size(tally%results, 2)
n = size(tally%results, 3)
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,:,:) = tally%results(:,:)%sum
tally_temp(2,:,:) = tally%results(:,:)%sum_sq
tally_temp(1,:,:) = tally%results(RESULT_SUM,:,:)
tally_temp(2,:,:) = tally%results(RESULT_SUM_SQ,:,:)
if (master) then
tally_group = open_group(tallies_group, "tally " // &
@ -592,20 +582,20 @@ contains
! At the end of the simulation, store the results back in the
! regular TallyResults array
if (current_batch == n_max_batches .or. satisfy_triggers) then
tally%results(:,:)%sum = tally_temp(1,:,:)
tally%results(:,:)%sum_sq = tally_temp(2,:,:)
tally%results(RESULT_SUM,:,:) = tally_temp(1,:,:)
tally%results(RESULT_SUM_SQ,:,:) = tally_temp(2,:,:)
end if
! Put in temporary tally result
allocate(tallyresult_temp(m,n))
tallyresult_temp(:,:)%sum = tally_temp(1,:,:)
tallyresult_temp(:,:)%sum_sq = tally_temp(2,:,:)
allocate(dummy_tally % results(3,m,n))
dummy_tally % results(RESULT_SUM,:,:) = tally_temp(1,:,:)
dummy_tally % results(RESULT_SUM_SQ,:,:) = tally_temp(2,:,:)
! Write reduced tally results to file
call write_dataset(tally_group, "results", tally%results)
call dummy_tally % write_results_hdf5(tally_group)
! Deallocate temporary tally result
deallocate(tallyresult_temp)
deallocate(dummy_tally % results)
else
! Receive buffer not significant at other processors
#ifdef MPI
@ -771,7 +761,7 @@ contains
call read_dataset(n_realizations, file_id, "n_realizations", indep=.true.)
! Read global tally data
call read_dataset(file_id, "global_tallies", global_tallies)
call read_dataset(global_tallies, file_id, "global_tallies")
! Check if tally results are present
tallies_group = open_group(file_id, "tallies")
@ -787,7 +777,7 @@ contains
! Read sum, sum_sq, and N for each bin
tally_group = open_group(tallies_group, "tally " // &
trim(to_str(tally % id)))
call read_dataset(tally_group, "results", tally % results)
call tally % read_results_hdf5(tally_group)
call read_dataset(tally % n_realizations, tally_group, &
"n_realizations")
call close_group(tally_group)

201
src/tally.F90
View file

@ -1,5 +1,7 @@
module tally
use, intrinsic :: ISO_C_BINDING
#ifdef MPI
use message_passing
#endif
@ -18,7 +20,6 @@ module tally
use output, only: header
use particle_header, only: LocalCoord, Particle
use string, only: to_str
use tally_header, only: TallyResult
use tally_filter
implicit none
@ -1965,8 +1966,8 @@ contains
score = score * calc_pn(t % moment_order(i), p % mu)
endif
!$omp atomic
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
t % results(RESULT_VALUE, score_index, filter_index) = &
t % results(RESULT_VALUE, score_index, filter_index) + score
case(SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN)
@ -1982,10 +1983,9 @@ contains
! multiply score by the angular flux moments and store
!$omp critical (score_general_scatt_yn)
t % results(score_index: score_index + num_nm - 1, filter_index) &
% value = t &
% results(score_index: score_index + num_nm - 1, filter_index)&
% value &
t % results(RESULT_VALUE, score_index: score_index + num_nm - 1, &
filter_index) = t % results(RESULT_VALUE, &
score_index: score_index + num_nm - 1, filter_index) &
+ score * calc_pn(n, p % mu) * calc_rn(n, p % last_uvw)
!$omp end critical (score_general_scatt_yn)
end do
@ -2011,10 +2011,9 @@ contains
! multiply score by the angular flux moments and store
!$omp critical (score_general_flux_tot_yn)
t % results(score_index: score_index + num_nm - 1, filter_index) &
% value = t &
% results(score_index: score_index + num_nm - 1, filter_index)&
% value &
t % results(RESULT_VALUE, score_index: score_index + num_nm - 1, &
filter_index) = t % results(RESULT_VALUE, &
score_index: score_index + num_nm - 1, filter_index) &
+ score * calc_rn(n, uvw)
!$omp end critical (score_general_flux_tot_yn)
end do
@ -2031,8 +2030,8 @@ contains
! get the score and tally it
!$omp atomic
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value &
t % results(RESULT_VALUE, score_index, filter_index) = &
t % results(RESULT_VALUE, score_index, filter_index) &
+ score * calc_pn(n, p % mu)
end do
i = i + t % moment_order(i)
@ -2040,8 +2039,8 @@ contains
case default
!$omp atomic
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
t % results(RESULT_VALUE, score_index, filter_index) = &
t % results(RESULT_VALUE, score_index, filter_index) + score
end select
@ -2458,8 +2457,8 @@ contains
! Add score to tally
!$omp atomic
t % results(i_score, i_filter) % value = &
t % results(i_score, i_filter) % value + score * filter_weight
t % results(RESULT_VALUE, i_score, i_filter) = &
t % results(RESULT_VALUE, i_score, i_filter) + score * filter_weight
! Case for tallying delayed emissions
else if (score_bin == SCORE_DELAYED_NU_FISSION .and. g /= 0) then
@ -2498,8 +2497,8 @@ contains
! Add score to tally
!$omp atomic
t % results(i_score, i_filter) % value = &
t % results(i_score, i_filter) % value + score * filter_weight
t % results(RESULT_VALUE, i_score, i_filter) = &
t % results(RESULT_VALUE, i_score, i_filter) + score * filter_weight
end if
end if
end do
@ -2536,8 +2535,8 @@ contains
filter_weight = product(filter_weights(:size(t % filters)))
!$omp atomic
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score * filter_weight
t % results(RESULT_VALUE, score_index, filter_index) = &
t % results(RESULT_VALUE, score_index, filter_index) + score * filter_weight
! reset original delayed group bin
matching_bins(t % find_filter(FILTER_DELAYEDGROUP)) = bin_original
@ -2997,8 +2996,8 @@ contains
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
!$omp atomic
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
t % results(RESULT_VALUE, 1, filter_index) = &
t % results(RESULT_VALUE, 1, filter_index) + p % wgt
end if
! Inward current on d1 min surface
@ -3033,8 +3032,8 @@ contains
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
!$omp atomic
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
t % results(RESULT_VALUE, 1, filter_index) = &
t % results(RESULT_VALUE, 1, filter_index) + p % wgt
ijk0(d1) = ijk0(d1) - 1
end if
@ -3053,8 +3052,8 @@ contains
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
!$omp atomic
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
t % results(RESULT_VALUE, 1, filter_index) = &
t % results(RESULT_VALUE, 1, filter_index) + p % wgt
end if
! Inward current on d1 max surface
@ -3089,8 +3088,8 @@ contains
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
!$omp atomic
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
t % results(RESULT_VALUE, 1, filter_index) = &
t % results(RESULT_VALUE, 1, filter_index) + p % wgt
ijk0(d1) = ijk0(d1) + 1
end if
@ -3116,9 +3115,10 @@ contains
subroutine synchronize_tallies()
integer :: i
real(8) :: k_col ! Copy of batch collision estimate of keff
real(8) :: k_abs ! Copy of batch absorption estimate of keff
real(8) :: k_tra ! Copy of batch tracklength estimate of keff
real(C_DOUBLE) :: k_col ! Copy of batch collision estimate of keff
real(C_DOUBLE) :: k_abs ! Copy of batch absorption estimate of keff
real(C_DOUBLE) :: k_tra ! Copy of batch tracklength estimate of keff
real(C_DOUBLE) :: val
#ifdef MPI
! Combine tally results onto master process
@ -3142,9 +3142,9 @@ contains
if (run_mode == MODE_EIGENVALUE) then
if (active_batches) then
! Accumulate products of different estimators of k
k_col = global_tallies(K_COLLISION) % value / total_weight
k_abs = global_tallies(K_ABSORPTION) % value / total_weight
k_tra = global_tallies(K_TRACKLENGTH) % value / total_weight
k_col = global_tallies(RESULT_VALUE, K_COLLISION) / total_weight
k_abs = global_tallies(RESULT_VALUE, K_ABSORPTION) / total_weight
k_tra = global_tallies(RESULT_VALUE, K_TRACKLENGTH) / total_weight
k_col_abs = k_col_abs + k_col * k_abs
k_col_tra = k_col_tra + k_col * k_tra
k_abs_tra = k_abs_tra + k_abs * k_tra
@ -3152,7 +3152,14 @@ contains
end if
! Accumulate results for global tallies
call accumulate_result(global_tallies)
do i = 1, size(global_tallies, 2)
val = global_tallies(RESULT_VALUE, i)/total_weight
global_tallies(RESULT_VALUE, i) = ZERO
global_tallies(RESULT_SUM, i) = global_tallies(RESULT_SUM, i) + val
global_tallies(RESULT_SUM_SQ, i) = &
global_tallies(RESULT_SUM_SQ, i) + val*val
end do
end if
end subroutine synchronize_tallies
@ -3182,7 +3189,7 @@ contains
allocate(tally_temp(m,n))
tally_temp = t % results(:,:) % value
tally_temp = t % results(RESULT_VALUE,:,:)
if (master) then
! The MPI_IN_PLACE specifier allows the master to copy values into
@ -3191,35 +3198,35 @@ contains
MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
! Transfer values to value on master
t % results(:,:) % value = tally_temp
t % results(RESULT_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)
! Reset value on other processors
t % results(:,:) % value = 0
t % results(RESULT_VALUE,:,:) = ZERO
end if
deallocate(tally_temp)
end do
! Copy global tallies into array to be reduced
global_temp = global_tallies(:) % value
global_temp = global_tallies(RESULT_VALUE, :)
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, global_temp, N_GLOBAL_TALLIES, &
MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
! Transfer values back to global_tallies on master
global_tallies(:) % value = global_temp
global_tallies(RESULT_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)
! Reset value on other processors
global_tallies(:) % value = ZERO
global_tallies(RESULT_VALUE, :) = ZERO
end if
! We also need to determine the total starting weight of particles from the
@ -3244,6 +3251,9 @@ contains
type(TallyObject), intent(inout) :: t
integer :: i, j
real(C_DOUBLE) :: val
! Increment number of realizations
if (reduce_tallies) then
t % n_realizations = t % n_realizations + 1
@ -3251,92 +3261,59 @@ contains
t % n_realizations = t % n_realizations + n_procs
end if
! Accumulate each TallyResult
call accumulate_result(t % results)
! Accumulate each result
do j = 1, size(t % results, 3)
do i = 1, size(t % results, 2)
val = t % results(RESULT_VALUE, i, j)/total_weight
t % results(RESULT_VALUE, i, j) = ZERO
t % results(RESULT_SUM, i, j) = &
t % results(RESULT_SUM, i, j) + val
t % results(RESULT_SUM_SQ, i, j) = &
t % results(RESULT_SUM_SQ, i, j) + val*val
end do
end do
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 TallyResults
! respectively
! tally and stores them in the RESULT_SUM and RESULT_SUM_SQ positions
!===============================================================================
subroutine tally_statistics()
integer :: i ! index in tallies array
type(TallyObject), pointer :: t
! Calculate statistics for user-defined tallies
do i = 1, n_tallies
t => tallies(i)
call statistics_result(t % results, t % n_realizations)
end do
! Calculate statistics for global tallies
call statistics_result(global_tallies, n_realizations)
end subroutine tally_statistics
!===============================================================================
! ACCUMULATE_RESULT accumulates results from many histories (or many generations)
! into a single realization of a random variable.
!===============================================================================
elemental subroutine accumulate_result(this)
type(TallyResult), intent(inout) :: this
real(8) :: val
! 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 = this % value/total_weight
this % sum = this % sum + val
this % sum_sq = this % sum_sq + val*val
! Reset the single batch estimate
this % value = ZERO
end subroutine accumulate_result
!===============================================================================
! STATISTICS_RESULT determines the sample mean and the standard deviation of the
! mean for a TallyResult.
!===============================================================================
elemental subroutine statistics_result(this, n)
type(TallyResult), intent(inout) :: this
integer, intent(in) :: n
integer :: j, k ! score/filter indices
integer :: n ! number of realizations
! 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.
this % sum = this % sum/n
this % sum_sq = sqrt((this % sum_sq/n - this % sum * &
this % sum) / (n - 1))
do i = 1, n_tallies
n = tallies(i) % n_realizations
end subroutine statistics_result
associate (r => tallies(i) % results)
do k = 1, size(r, 3)
do j = 1, size(r, 2)
r(RESULT_SUM, j, k) = r(RESULT_SUM, j, k) / n
r(RESULT_SUM_SQ, j, k) = sqrt((r(RESULT_SUM_SQ, j, k)/n - &
r(RESULT_SUM, j, k) * r(RESULT_SUM, j, k))/(n - 1))
end do
end do
end associate
end do
!===============================================================================
! RESET_RESULT zeroes out the value and accumulated sum and sum-squared for a
! single TallyResult.
!===============================================================================
elemental subroutine reset_result(this)
type(TallyResult), intent(inout) :: this
this % value = ZERO
this % sum = ZERO
this % sum_sq = ZERO
end subroutine reset_result
! Calculate statistics for global tallies
n = n_realizations
associate (r => global_tallies)
do j = 1, size(r, 2)
r(RESULT_SUM, j) = r(RESULT_SUM, j) / n
r(RESULT_SUM_SQ, j) = sqrt((r(RESULT_SUM_SQ, j)/n - &
r(RESULT_SUM, j) * r(RESULT_SUM, j))/(n - 1))
end do
end associate
end subroutine tally_statistics
!===============================================================================
! SETUP_ACTIVE_USERTALLIES

91
src/tally_header.F90
View file

@ -1,23 +1,15 @@
module tally_header
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants, only: NONE, N_FILTER_TYPES
use tally_filter_header, only: TallyFilterContainer
use trigger_header, only: TriggerObject
use, intrinsic :: ISO_C_BINDING
implicit none
!===============================================================================
! TALLYRESULT provides accumulation of results in a particular tally bin
!===============================================================================
type, bind(C) :: TallyResult
real(C_DOUBLE) :: value = 0.
real(C_DOUBLE) :: sum = 0.
real(C_DOUBLE) :: sum_sq = 0.
end type TallyResult
!===============================================================================
! TALLYOBJECT describes a user-specified tally. The region of phase space to
! tally in is given by the TallyFilters and the results are stored in a
@ -68,7 +60,7 @@ module tally_header
integer :: total_filter_bins
integer :: total_score_bins
type(TallyResult), allocatable :: results(:,:)
real(C_DOUBLE), allocatable :: results(:,:,:)
! reset property - allows a tally to be reset after every batch
logical :: reset = .false.
@ -79,6 +71,79 @@ module tally_header
! Tally precision triggers
integer :: n_triggers = 0 ! # of triggers
type(TriggerObject), allocatable :: triggers(:) ! Array of triggers
contains
procedure :: write_results_hdf5
procedure :: read_results_hdf5
end type TallyObject
contains
subroutine write_results_hdf5(this, group_id)
class(TallyObject), intent(in) :: this
integer(HID_T), intent(in) :: group_id
integer :: hdf5_err
integer(HID_T) :: dset, dspace
integer(HID_T) :: memspace
integer(HSIZE_T) :: dims(3)
integer(HSIZE_T) :: dims_slab(3)
integer(HSIZE_T) :: offset(3) = [1,0,0]
! Create file dataspace
dims_slab(:) = shape(this % results)
dims_slab(1) = 2
call h5screate_simple_f(3, dims_slab, dspace, hdf5_err)
! Create memory dataspace that contains only SUM and SUM_SQ values
dims(:) = shape(this % results)
call h5screate_simple_f(3, dims, memspace, hdf5_err)
call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offset, dims_slab, &
hdf5_err)
! Create and write to dataset
call h5dcreate_f(group_id, "results", H5T_NATIVE_DOUBLE, dspace, dset, &
hdf5_err)
call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, this % results, dims_slab, &
hdf5_err, mem_space_id=memspace)
! Close identifiers
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(memspace, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
end subroutine write_results_hdf5
subroutine read_results_hdf5(this, group_id)
class(TallyObject), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
integer :: hdf5_err
integer(HID_T) :: dset, dspace
integer(HID_T) :: memspace
integer(HSIZE_T) :: dims(3)
integer(HSIZE_T) :: dims_slab(3)
integer(HSIZE_T) :: offset(3) = [1,0,0]
! Create file dataspace
dims_slab(:) = shape(this % results)
dims_slab(1) = 2
call h5screate_simple_f(3, dims_slab, dspace, hdf5_err)
! Create memory dataspace that contains only SUM and SUM_SQ values
dims(:) = shape(this % results)
call h5screate_simple_f(3, dims, memspace, hdf5_err)
call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, offset, dims_slab, &
hdf5_err)
! Create and write to dataset
call h5dopen_f(group_id, "results", dset, hdf5_err)
call h5dread_f(dset, H5T_NATIVE_DOUBLE, this % results, dims_slab, &
hdf5_err, mem_space_id=memspace)
! Close identifiers
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(memspace, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
end subroutine read_results_hdf5
end module tally_header

6
src/tally_initialize.F90
View file

@ -20,7 +20,8 @@ contains
subroutine configure_tallies()
! Allocate global tallies
allocate(global_tallies(N_GLOBAL_TALLIES))
allocate(global_tallies(3, N_GLOBAL_TALLIES))
global_tallies(:,:) = ZERO
call setup_tally_arrays()
@ -62,7 +63,8 @@ contains
t % total_score_bins = t % n_score_bins * t % n_nuclide_bins
! Allocate results array
allocate(t % results(t % total_score_bins, t % total_filter_bins))
allocate(t % results(3, t % total_score_bins, t % total_filter_bins))
t % results(:,:,:) = ZERO
end do TALLY_LOOP

16
src/trigger.F90
View file

@ -432,17 +432,19 @@ contains
real(8), intent(inout) :: rel_err ! tally relative error
integer, intent(in) :: score_index ! tally results score index
integer, intent(in) :: filter_index ! tally results filter index
integer :: n ! number of realizations
real(8) :: mean ! tally mean
type(TallyResult) :: tally_result ! pointer to TallyResult
type(TallyObject), pointer :: t ! tally pointer
type(TallyObject), intent(in) :: t ! tally
integer :: n ! number of realizations
real(8) :: mean ! tally mean
real(8) :: tally_sum, tally_sum_sq ! results for a single tally bin
n = t % n_realizations
tally_result = t % results(score_index, filter_index)
tally_sum = t % results(RESULT_SUM, score_index, filter_index)
tally_sum_sq = t % results(RESULT_SUM_SQ, score_index, filter_index)
! Compute the tally mean and standard deviation
mean = tally_result % sum / n
std_dev = sqrt((tally_result % sum_sq / n - mean * mean) / (n - 1))
mean = tally_sum / n
std_dev = sqrt((tally_sum_sq / n - mean * mean) / (n - 1))
! Compute the relative error if the mean is non-zero
if (mean == ZERO) then