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

Merge pull request #8 from paulromano/python-api-pr

Browse source 
A few updates for Python API
This commit is contained in:
Will Boyd 2015-04-17 10:36:54 -04:00
commit 83db4c5650
44 changed files with 1291 additions and 1499 deletions
Download patch file Download diff file Expand all files Collapse all files

3
.gitignore vendored
View file

@ -11,6 +11,9 @@
# OpenMC executable
src/openmc
# Inputs generated from Python API
examples/python/**/*.xml
# emacs backups
*~

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

@ -174,11 +174,15 @@ should be performed. It has the following attributes/sub-elements:
-------------------------
The ``<energy_grid>`` element determines the treatment of the energy grid during
a simulation. The valid options are "nuclide" and "logarithm". Setting this
element to "nuclide" will cause OpenMC to use a nuclide's energy grid when
determining what points to interpolate between for determining cross sections
(i.e. non-unionized energy grid). Setting this element to "logarithm" causes
OpenMC to use a logarithmic mapping technique described in LA-UR-14-24530_.
a simulation. The valid options are "nuclide", "logarithm", and
"material-union". Setting this element to "nuclide" will cause OpenMC to use a
nuclide's energy grid when determining what points to interpolate between for
determining cross sections (i.e. non-unionized energy grid). Setting this
element to "logarithm" causes OpenMC to use a logarithmic mapping technique
described in LA-UR-14-24530_. Setting this element to "material-union" will
cause OpenMC to create energy grids that are unionized material-by-material and
use these grids when determining the energy-cross section pairs to interpolate
cross section values between.
*Default*: logarithm

2
src/CMakeLists.txt
View file

@ -37,7 +37,7 @@ endif()
set(MPI_ENABLED FALSE)
set(HDF5_ENABLED FALSE)
if($ENV{FC} MATCHES "mpi.*")
if($ENV{FC} MATCHES "mpi[^/]*$")
message("-- Detected MPI wrapper: $ENV{FC}")
add_definitions(-DMPI)
set(MPI_ENABLED TRUE)

5
src/ace_header.F90
View file

@ -103,7 +103,7 @@ module ace_header
! Energy grid information
integer :: n_grid ! # of nuclide grid points
integer, allocatable :: grid_index(:) ! union grid pointers / log grid mapping
integer, allocatable :: grid_index(:) ! log grid mapping indices
real(8), allocatable :: energy(:) ! energy values corresponding to xs
! Microscopic cross sections
@ -372,9 +372,6 @@ module ace_header
integer :: i ! Loop counter
if (allocated(this % grid_index)) &
deallocate(this % grid_index)
if (allocated(this % energy)) &
deallocate(this % energy, this % total, this % elastic, &
& this % fission, this % nu_fission, this % absorption)

5
src/constants.F90
View file

@ -368,8 +368,9 @@ module constants
! Energy grid methods
integer, parameter :: &
GRID_NUCLIDE = 1, & ! non-unionized energy grid
GRID_LOGARITHM = 2 ! logarithmic mapping
GRID_NUCLIDE = 1, & ! unique energy grid for each nuclide
GRID_MAT_UNION = 2, & ! material union grids with pointers
GRID_LOGARITHM = 3 ! lethargy mapping
! Running modes
integer, parameter :: &

75
src/cross_section.F90
View file

@ -15,6 +15,9 @@ module cross_section
implicit none
save
integer :: union_grid_index
!$omp threadprivate(union_grid_index)
contains
!===============================================================================
@ -36,11 +39,11 @@ contains
!$omp threadprivate(mat)
! Set all material macroscopic cross sections to zero
material_xs % total = ZERO
material_xs % elastic = ZERO
material_xs % absorption = ZERO
material_xs % fission = ZERO
material_xs % nu_fission = ZERO
material_xs % total = ZERO
material_xs % elastic = ZERO
material_xs % absorption = ZERO
material_xs % fission = ZERO
material_xs % nu_fission = ZERO
material_xs % kappa_fission = ZERO
! Exit subroutine if material is void
@ -48,6 +51,10 @@ contains
mat => materials(p % material)
! Find energy index on global or material unionized grid
if (grid_method == GRID_MAT_UNION) &
call find_energy_index(p % E, p % material)
! Determine if this material has S(a,b) tables
check_sab = (mat % n_sab > 0)
@ -88,9 +95,9 @@ contains
! Calculate microscopic cross section for this nuclide
if (p % E /= micro_xs(i_nuclide) % last_E) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E)
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i)
else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E)
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i)
end if
! ========================================================================
@ -131,24 +138,32 @@ contains
! given index in the nuclides array at the energy of the given particle
!===============================================================================
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E)
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
integer :: i_grid ! index on nuclide energy grid
integer :: i_low, i_high ! bounding indices from logarithmic mapping
integer :: u ! index into logarithmic mapping array
integer, intent(in) :: i_mat ! index into materials array
integer, intent(in) :: i_nuc_mat ! index into nuclides array for a material
integer :: i_grid ! index on nuclide energy grid
integer :: i_low ! lower logarithmic mapping index
integer :: i_high ! upper logarithmic mapping index
integer :: u ! index into logarithmic mapping array
real(8), intent(in) :: E ! energy
real(8) :: f ! interp factor on nuclide energy grid
type(Nuclide), pointer, save :: nuc => null()
!$omp threadprivate(nuc)
type(Nuclide), pointer, save :: nuc => null()
type(Material), pointer, save :: mat => null()
!$omp threadprivate(nuc, mat)
! Set pointer to nuclide
! Set pointer to nuclide and material
nuc => nuclides(i_nuclide)
mat => materials(i_mat)
! Determine index on nuclide energy grid
select case (grid_method)
case (GRID_MAT_UNION)
i_grid = mat % nuclide_grid_index(i_nuc_mat, union_grid_index)
case (GRID_LOGARITHM)
! Determine the energy grid index using a logarithmic mapping to reduce
! the energy range over which a binary search needs to be performed
@ -173,7 +188,7 @@ contains
! Perform binary search on the nuclide energy grid in order to determine
! which points to interpolate between
if (E < nuc % energy(1)) then
if (E <= nuc % energy(1)) then
i_grid = 1
elseif (E > nuc % energy(nuc % n_grid)) then
i_grid = nuc % n_grid - 1
@ -506,6 +521,32 @@ contains
end subroutine calculate_urr_xs
!===============================================================================
! FIND_ENERGY_INDEX determines the index on the union energy grid at a certain
! energy
!===============================================================================
subroutine find_energy_index(E, i_mat)
real(8), intent(in) :: E ! energy of particle
integer, intent(in) :: i_mat ! material index
type(Material), pointer, save :: mat => null() ! pointer to current material
!$omp threadprivate(mat)
mat => materials(i_mat)
! if the energy is outside of energy grid range, set to first or last
! index. Otherwise, do a binary search through the union energy grid.
if (E <= mat % e_grid(1)) then
union_grid_index = 1
elseif (E > mat % e_grid(mat % n_grid)) then
union_grid_index = mat % n_grid - 1
else
union_grid_index = binary_search(mat % e_grid, mat % n_grid, E)
end if
end subroutine find_energy_index
!===============================================================================
! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section
! for a given nuclide at the trial relative energy used in resonance scattering

155
src/energy_grid.F90
View file

@ -1,6 +1,9 @@
module energy_grid
use global
use list_header, only: ListReal
use output, only: write_message
implicit none
@ -10,6 +13,53 @@ module energy_grid
contains
!===============================================================================
! UNIONIZED_GRID creates a unionized energy grid, for the entire problem or for
! each material, composed of the grids from each nuclide in the entire problem,
! or each material, respectively. Right now, the grid for each nuclide is added
! into a linked list one at a time with an effective insertion sort. Could be
! done with a hash for all energy points and then a quicksort at the end (what
! hash function to use?)
!===============================================================================
subroutine unionized_grid()
integer :: i ! index in nuclides array
integer :: j ! index in materials array
type(ListReal), pointer, save :: list => null()
type(Nuclide), pointer, save :: nuc => null()
type(Material), pointer, save :: mat => null()
!$omp threadprivate(list, nuc, mat)
call write_message("Creating unionized energy grid...", 5)
! add grid points for each nuclide in the material
do j = 1, n_materials
mat => materials(j)
do i = 1, mat % n_nuclides
nuc => nuclides(mat % nuclide(i))
call add_grid_points(list, nuc % energy)
end do
! set size of unionized material energy grid
mat % n_grid = list % size()
! create allocated array from linked list
allocate(mat % e_grid(mat % n_grid))
do i = 1, mat % n_grid
mat % e_grid(i) = list % get_item(i)
end do
! delete linked list and dictionary
call list % clear()
deallocate(list)
end do
! Set pointers to unionized energy grid for each nuclide
call grid_pointers()
end subroutine unionized_grid
!===============================================================================
! LOGARITHMIC_GRID determines a logarithmic mapping for energies to bounding
! indices on a nuclide energy grid
@ -59,4 +109,109 @@ contains
end subroutine logarithmic_grid
!===============================================================================
! ADD_GRID_POINTS adds energy points from the 'energy' array into a linked list
! of points already stored from previous arrays.
!===============================================================================
subroutine add_grid_points(list, energy)
type(ListReal), pointer :: list
real(8), intent(in) :: energy(:)
integer :: i ! index in energy array
integer :: n ! size of energy array
integer :: current ! current index
real(8) :: E ! actual energy value
i = 1
n = size(energy)
! If the original list is empty, we need to allocate the first element and
! store first energy point
if (.not. associated(list)) then
allocate(list)
do i = 1, n
call list % append(energy(i))
end do
return
end if
! Set current index to beginning of the list
current = 1
do while (i <= n)
E = energy(i)
! If we've reached the end of the grid energy list, add the remaining
! energy points to the end
if (current > list % size()) then
! Finish remaining energies
do while (i <= n)
call list % append(energy(i))
i = i + 1
end do
exit
end if
if (E < list % get_item(current)) then
! Insert new energy in this position
call list % insert(current, E)
! Advance index in linked list and in new energy grid
i = i + 1
current = current + 1
elseif (E == list % get_item(current)) then
! Found the exact same energy, no need to store duplicates so just
! skip and move to next index
i = i + 1
current = current + 1
else
current = current + 1
end if
end do
end subroutine add_grid_points
!===============================================================================
! GRID_POINTERS creates an array of pointers (ints) for each nuclide to link
! each point on the nuclide energy grid to one on a unionized energy grid
!===============================================================================
subroutine grid_pointers()
integer :: i ! loop index for nuclides
integer :: j ! loop index for nuclide energy grid
integer :: k ! loop index for materials
integer :: index_e ! index on union energy grid
real(8) :: union_energy ! energy on union grid
real(8) :: energy ! energy on nuclide grid
type(Nuclide), pointer :: nuc => null()
type(Material), pointer :: mat => null()
do k = 1, n_materials
mat => materials(k)
allocate(mat % nuclide_grid_index(mat % n_nuclides, mat % n_grid))
do i = 1, mat % n_nuclides
nuc => nuclides(mat % nuclide(i))
index_e = 1
energy = nuc % energy(index_e)
do j = 1, mat % n_grid
union_energy = mat % e_grid(j)
if (union_energy >= energy .and. index_e < nuc % n_grid) then
index_e = index_e + 1
energy = nuc % energy(index_e)
end if
mat % nuclide_grid_index(i,j) = index_e - 1
end do
end do
end do
end subroutine grid_pointers
end module energy_grid

1
src/global.F90
View file

@ -218,6 +218,7 @@ module global
type(Timer) :: time_total ! timer for total run
type(Timer) :: time_initialize ! timer for initialization
type(Timer) :: time_read_xs ! timer for reading cross sections
type(Timer) :: time_unionize ! timer for material xs-energy grid union
type(Timer) :: time_bank ! timer for fission bank synchronization
type(Timer) :: time_bank_sample ! timer for fission bank sampling
type(Timer) :: time_bank_sendrecv ! timer for fission bank SEND/RECV

71
src/hdf5_interface.F90
View file

@ -68,7 +68,6 @@ module hdf5_interface
module procedure hdf5_read_integer_4Darray
module procedure hdf5_read_long
module procedure hdf5_read_string
module procedure hdf5_read_string_1Darray
#ifdef MPI
module procedure hdf5_read_double_parallel
module procedure hdf5_read_double_1Darray_parallel
@ -82,7 +81,6 @@ module hdf5_interface
module procedure hdf5_read_integer_4Darray_parallel
module procedure hdf5_read_long_parallel
module procedure hdf5_read_string_parallel
! module procedure hdf5_read_string_1Darray_parallel
#endif
end interface hdf5_read_data
@ -208,7 +206,7 @@ contains
logical :: status ! does the group exist
! Check if group exists
call h5ltpath_valid_f(hdf5_fh, trim(group), .true., status, hdf5_err)
call h5ltpath_valid_f(hdf5_fh, trim(group), .true., status, hdf5_err)
! Either create or open group
if (status) then
@ -259,7 +257,7 @@ contains
integer(HID_T), intent(in) :: group ! name of group
character(*), intent(in) :: name ! name of data
integer, intent(inout) :: buffer ! read data to here
integer, intent(inout) :: buffer ! read data to here
integer :: buffer_copy(1) ! need an array for read
@ -463,7 +461,7 @@ contains
integer(HID_T), intent(in) :: group ! name of group
character(*), intent(in) :: name ! name of data
real(8), intent(inout) :: buffer ! read data to here
real(8), intent(inout) :: buffer ! read data to here
real(8) :: buffer_copy(1) ! need an array for read
@ -728,7 +726,7 @@ contains
call h5dcreate_f(group, name, H5T_STRING, dspace, dset, hdf5_err)
! Set up dimesnions of string to write
dims2 = (/length, 1/) ! full array of strings to write
dims2 = (/length, 1/) ! full array of strings to write
dims1(1) = length ! length of string
! Copy over string buffer to a rank 1 array
@ -762,7 +760,7 @@ contains
! character(len=length, kind=c_char), pointer :: chr_ptr
! f_ptr = c_loc(buf_ptr(1))
! call h5dread_f(dset, H5T_STRING, f_ptr, hdf5_err, xfer_prp=plist)
! call c_f_pointer(buf_ptr(1), chr_ptr)
! call c_f_pointer(buf_ptr(1), chr_ptr)
! buffer = chr_ptr
! nullify(chr_ptr)
@ -788,39 +786,6 @@ contains
end subroutine hdf5_read_string
!===============================================================================
! HDF5_READ_STRING_1DARRAY reads string data
!===============================================================================
subroutine hdf5_read_string_1Darray(group, name, buffer, length, rank)
integer(HID_T), intent(in) :: group ! name of group
character(*), intent(in) :: name ! name of data
integer, intent(in) :: length ! length of strings
integer, intent(in) :: rank ! number of strings
character(*), intent(inout) :: buffer(length,rank) ! read data to here
! Open dataset
call h5dopen_f(group, name, dset, hdf5_err)
! Get dataspace to read
call h5dget_space_f(dset, dspace, hdf5_err)
! Set dimensions
dims2 = (/length, rank/)
dims1(1) = length
! Read in the data
!call h5ltread_dataset_string_f(dset, H5T_STRING, buffer, dims2, dims1, hdf5_err, &
! mem_space_id=dspace, xfer_prp = plist)
!call h5ltread_dataset_f(group, name, dims, dims2, H5T_STRING, buffer, hdf5_err)
! Close dataset
call h5dclose_f(dset, hdf5_err)
end subroutine hdf5_read_string_1Darray
!===============================================================================
! HDF5_WRITE_ATTRIBUTE_STRING writes a string attribute to a variables
!===============================================================================
@ -873,7 +838,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -951,7 +916,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1031,7 +996,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1112,7 +1077,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1194,7 +1159,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1273,7 +1238,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1351,7 +1316,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1431,7 +1396,7 @@ contains
f_ptr = c_loc(buffer(1,1))
call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1512,7 +1477,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1594,7 +1559,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1674,7 +1639,7 @@ contains
f_ptr = c_loc(buffer)
call h5dwrite_f(dset, long_type, f_ptr, hdf5_err, xfer_prp=plist)
! Close all
! Close all
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5pclose_f(plist, hdf5_err)
@ -1761,7 +1726,7 @@ contains
call h5dcreate_f(group, name, H5T_STRING, dspace, dset, hdf5_err)
! Set up dimesnions of string to write
dims2 = (/length, 1/) ! full array of strings to write
dims2 = (/length, 1/) ! full array of strings to write
dims1(1) = length ! length of string
! Copy over string buffer to a rank 1 array
@ -1797,7 +1762,7 @@ contains
! character(len=length, kind=c_char), pointer :: chr_ptr
! f_ptr = c_loc(buf_ptr(1))
! call h5dread_f(dset, H5T_STRING, f_ptr, hdf5_err, xfer_prp=plist)
! call c_f_pointer(buf_ptr(1), chr_ptr)
! call c_f_pointer(buf_ptr(1), chr_ptr)
! buffer = chr_ptr
! nullify(chr_ptr)

15
src/initialize.F90
View file

@ -4,7 +4,7 @@ module initialize
use bank_header, only: Bank
use constants
use dict_header, only: DictIntInt, ElemKeyValueII
use energy_grid, only: logarithmic_grid, grid_method
use energy_grid, only: logarithmic_grid, grid_method, unionized_grid
use error, only: fatal_error, warning
use geometry, only: neighbor_lists
use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,&
@ -109,10 +109,17 @@ contains
! Create linked lists for multiple instances of the same nuclide
call same_nuclide_list()
! Construct logarithmic energy grid for cross-sections
if (grid_method == GRID_LOGARITHM) then
! Construct unionized or log energy grid for cross-sections
select case (grid_method)
case (GRID_NUCLIDE)
continue
case (GRID_MAT_UNION)
call time_unionize % start()
call unionized_grid()
call time_unionize % stop()
case (GRID_LOGARITHM)
call logarithmic_grid()
end if
end select
! Allocate and setup tally stride, matching_bins, and tally maps
call configure_tallies()

7
src/input_xml.F90
View file

@ -213,8 +213,11 @@ contains
select case (trim(temp_str))
case ('nuclide')
grid_method = GRID_NUCLIDE
case ('union')
call fatal_error("Union energy grid is no longer supported.")
case ('material-union', 'union')
grid_method = GRID_MAT_UNION
if (trim(temp_str) == 'union') &
call warning('Energy grids will be unionized by material. Global&
& energy grid unionization is no longer an allowed option.')
case ('logarithm', 'logarithmic', 'log')
grid_method = GRID_LOGARITHM
case default

7
src/material_header.F90
View file

@ -13,6 +13,13 @@ module material_header
real(8) :: density ! total atom density in atom/b-cm
real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm
! Energy grid information
integer :: n_grid ! # of union material grid points
real(8), allocatable :: e_grid(:) ! union material grid energies
! Unionized energy grid information
integer, allocatable :: nuclide_grid_index(:,:) ! nuclide e_grid pointers
! S(a,b) data references
integer :: n_sab = 0 ! number of S(a,b) tables
integer, allocatable :: i_sab_nuclides(:) ! index of corresponding nuclide

24
src/output.F90
View file

@ -1612,22 +1612,26 @@ contains
! write global tallies
if (n_realizations > 1) 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
if (n_realizations > 3) write(ou,102) "Combined k-effective", k_combined
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
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
else
if (master) call warning("Could not compute uncertainties -- only one &
&active batch simulated!")
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
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
end if
write(ou,103) "Leakage Fraction", global_tallies(LEAKAGE) % sum
end if
write(ou,*)

2
src/relaxng/settings.rnc
View file

@ -27,7 +27,7 @@ element settings {
(element weight_avg { xsd:double } | attribute weight_avg { xsd:double })?
}? &
element energy_grid { ( "nuclide" | "log" | "logarithm" | "logarithmic" ) }? &
element energy_grid { ( "nuclide" | "log" | "logarithm" | "logarithmic" | "material-union" | "union" ) }? &
element entropy {
(element dimension { list { xsd:int+ } } |

2
src/relaxng/settings.rng
View file

@ -106,6 +106,8 @@
<value>log</value>
<value>logarithm</value>
<value>logarithmic</value>
<value>material-union</value>
<value>union</value>
</choice>
</element>
</optional>

34
src/state_point.F90
View file

@ -39,17 +39,13 @@ contains
subroutine write_state_point()
character(MAX_FILE_LEN) :: filename
integer :: i, j, k, m
integer :: n_x, n_y, n_z
character(11), allocatable :: name_array(:)
integer :: i, j, k
integer, allocatable :: id_array(:)
integer, allocatable :: key_array(:)
type(StructuredMesh), pointer :: mesh => null()
type(TallyObject), pointer :: tally => null()
type(ElemKeyValueII), pointer :: current => null()
type(ElemKeyValueII), pointer :: next => null()
type(ElemKeyValueCI), pointer :: cur_nuclide => null()
type(ElemKeyValueCI), pointer :: next_nuclide => null()
type(StructuredMesh), pointer :: mesh
type(TallyObject), pointer :: tally
type(ElemKeyValueII), pointer :: current
type(ElemKeyValueII), pointer :: next
character(8) :: moment_name ! name of moment (e.g, P3)
integer :: n_order ! loop index for moment orders
integer :: nm_order ! loop index for Ynm moment orders
@ -306,14 +302,14 @@ contains
case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N)
moment_name = 'P' // to_str(tally % moment_order(k))
call sp % write_data(moment_name, "order" // trim(to_str(k)), &
group="tallies/tally " // trim(to_str(tally % id)) // &
group="tallies/tally " // trim(to_str(tally % id)) // &
"/moments")
k = k + 1
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
do n_order = 0, tally % moment_order(k)
moment_name = 'P' // trim(to_str(n_order))
call sp % write_data(moment_name, "order" // trim(to_str(k)), &
group="tallies/tally " // trim(to_str(tally % id)) // &
group="tallies/tally " // trim(to_str(tally % id)) // &
"/moments")
k = k + 1
end do
@ -325,7 +321,7 @@ contains
trim(to_str(nm_order))
call sp % write_data(moment_name, "order" // &
trim(to_str(k)), &
group="tallies/tally " // trim(to_str(tally % id)) // &
group="tallies/tally " // trim(to_str(tally % id)) // &
"/moments")
k = k + 1
end do
@ -504,9 +500,9 @@ contains
real(8) :: dummy ! temporary receive buffer for non-root reduces
#endif
integer, allocatable :: id_array(:)
type(ElemKeyValueII), pointer :: current => null()
type(ElemKeyValueII), pointer :: next => null()
type(TallyObject), pointer :: tally => null()
type(ElemKeyValueII), pointer :: current
type(ElemKeyValueII), pointer :: next
type(TallyObject), pointer :: tally
type(TallyResult), allocatable :: tallyresult_temp(:,:)
! ==========================================================================
@ -661,14 +657,10 @@ contains
integer, allocatable :: key_array(:)
integer :: curr_key
integer, allocatable :: temp_array(:)
integer, allocatable :: temp_array3D(:,:,:)
integer, allocatable :: temp_array4D(:,:,:,:)
logical :: source_present
real(8) :: l
real(8) :: real_array(3)
real(8), allocatable :: temp_real_array(:)
type(StructuredMesh), pointer :: mesh => null()
type(TallyObject), pointer :: tally => null()
type(StructuredMesh), pointer :: mesh
type(TallyObject), pointer :: tally
integer :: n_order ! loop index for moment orders
integer :: nm_order ! loop index for Ynm moment orders
character(8) :: moment_name ! name of moment (e.g, P3, Y-1,1)

1696
src/tally.F90
View file

File diff suppressed because it is too large Load diff

2
src/utils/openmc/__init__.py
View file

@ -1,5 +1,3 @@
#!/usr/bin/env python
from openmc.element import *
from openmc.geometry import *
from openmc.nuclide import *

3
src/utils/openmc/checkvalue.py
View file

@ -1,5 +1,3 @@
#!/usr/bin/env python
import numpy as np
@ -13,4 +11,3 @@ def is_float(val):
def is_string(val):
return isinstance(val, (str, np.str))

14
src/utils/openmc/clean_xml.py
View file

@ -1,12 +1,10 @@
#!/usr/bin/env python
def sort_xml_elements(tree):
# Retrieve all children of the root XML node in the tree
elements = tree.getchildren()
# Initialize empty lists for the sorted and comment elements
sorted_elements = list()
sorted_elements = []
# Initialize an empty set of tags (e.g., Surface, Cell, and Lattice)
tags = set()
@ -16,7 +14,7 @@ def sort_xml_elements(tree):
tags.add(element.tag)
# Initialize an empty list for the comment elements
comment_elements = list()
comment_elements = []
# Find the comment elements and record their ordering within the
# tree using a precedence with respect to the subsequent nodes
@ -40,7 +38,7 @@ def sort_xml_elements(tree):
continue
# Initialize an empty list of tuples to sort (id, element)
tagged_data = list()
tagged_data = []
# Retrieve the IDs for each of the elements
for element in tagged_elements:
@ -68,11 +66,11 @@ def sort_xml_elements(tree):
def clean_xml_indentation(element, level=0):
'''
"""
copy and paste from http://effbot.org/zone/elementent-lib.htm#prettyprint
it basically walks your tree and adds spaces and newlines so the tree is
printed in a nice way
'''
"""
i = "\n" + level*" "
@ -92,4 +90,4 @@ def clean_xml_indentation(element, level=0):
else:
if level and (not element.tail or not element.tail.strip()):
element.tail = i
element.tail = i

8
src/utils/openmc/cmfd.py
View file

@ -1,9 +1,9 @@
#!/usr/bin/env python
from xml.etree import ElementTree as ET
import numpy as np
from openmc.checkvalue import *
from openmc.clean_xml import *
from xml.etree import ElementTree as ET
import numpy as np
class CMFDMesh(object):
@ -583,4 +583,4 @@ class CMFDFile(object):
# Write the XML Tree to the cmfd.xml file
tree = ET.ElementTree(self._cmfd_file)
tree.write("cmfd.xml", xml_declaration=True,
encoding='utf-8', method="xml")
encoding='utf-8', method="xml")

2
src/utils/openmc/constants.py
View file

@ -1,5 +1,3 @@
#!/usr/bin/env python
"""Dictionaries of integer-to-string mappings from openmc/src/constants.F90"""
SURFACE_TYPES = {1: 'x-plane',

4
src/utils/openmc/element.py
View file

@ -1,5 +1,3 @@
#!/usr/bin/env python
from openmc.checkvalue import *
class Element(object):
@ -36,7 +34,7 @@ class Element(object):
def __hash__(self):
hashable = list()
hashable = []
hashable.append(self._name)
hashable.append(self._xs)
return hash(tuple(hashable))

14
src/utils/openmc/executor.py
View file

@ -1,11 +1,7 @@
#!/usr/bin/env python
from openmc.checkvalue import *
import subprocess
import os
FNULL = open(os.devnull, 'w')
from openmc.checkvalue import *
class Executor(object):
@ -36,7 +32,8 @@ class Executor(object):
subprocess.check_call('openmc -p', shell=True,
cwd=self._working_directory)
else:
subprocess.check_call('openmc -p', shell=True, stdout=FNULL,
subprocess.check_call('openmc -p', shell=True,
stdout=open(os.devnull, 'w'),
cwd=self._working_directory)
@ -71,5 +68,6 @@ class Executor(object):
subprocess.check_call(command, shell=True,
cwd=self._working_directory)
else:
subprocess.check_call(command, shell=True, stdout=FNULL,
cwd=self._working_directory)
subprocess.check_call(command, shell=True,
stdout=open(os.devnull, 'w'),
cwd=self._working_directory)

10
src/utils/openmc/geometry.py
View file

@ -1,8 +1,8 @@
#!/usr/bin/env python
from xml.etree import ElementTree as ET
import openmc
from openmc.clean_xml import *
from xml.etree import ElementTree as ET
def reset_auto_ids():
openmc.reset_auto_material_id()
@ -17,7 +17,7 @@ class Geometry(object):
# Initialize Geometry class attributes
self._root_universe = None
self._offsets = dict()
self._offsets = {}
def get_offset(self, path, filter_offset):
@ -63,7 +63,7 @@ class Geometry(object):
def get_all_nuclides(self):
nuclides = dict()
nuclides = {}
materials = self.get_all_materials()
for material in materials:
@ -159,4 +159,4 @@ class GeometryFile(object):
# Write the XML Tree to the materials.xml file
tree = ET.ElementTree(self._geometry_file)
tree.write("geometry.xml", xml_declaration=True,
encoding='utf-8', method="xml")
encoding='utf-8', method="xml")

38
src/utils/openmc/material.py
View file

@ -1,18 +1,17 @@
#!/usr/bin/env python
from collections import MappingView
from copy import deepcopy
import warnings
from xml.etree import ElementTree as ET
import numpy as np
import openmc
from openmc.checkvalue import *
from openmc.clean_xml import *
from xml.etree import ElementTree as ET
from collections import MappingView
from copy import deepcopy
import numpy as np
# A list of all IDs for all Materials created
MATERIAL_IDS = list()
MATERIAL_IDS = []
# A static variable for auto-generated Material IDs
AUTO_MATERIAL_ID = 10000
@ -20,7 +19,7 @@ AUTO_MATERIAL_ID = 10000
def reset_auto_material_id():
global AUTO_MATERIAL_ID, MATERIAL_IDS
AUTO_MATERIAL_ID = 10000
MATERIAL_IDS = list()
MATERIAL_IDS = []
# Units for density supported by OpenMC
@ -49,15 +48,15 @@ class Material(object):
# A dictionary of Nuclides
# Keys - Nuclide names
# Values - tuple (nuclide, percent, percent type)
self._nuclides = dict()
self._nuclides = {}
# A dictionary of Elements
# Keys - Element names
# Values - tuple (element, percent, percent type)
self._elements = dict()
self._elements = {}
# If specified, a list of tuples of (table name, xs identifier)
self._sab = list()
self._sab = []
# If true, the material will be initialized as distributed
self._convert_to_distrib_comps = False
@ -72,14 +71,13 @@ class Material(object):
def set_id(self, material_id=None):
global MATERIAL_IDS
global AUTO_MATERIAL_ID, MATERIAL_IDS
# If the Material already has an ID, remove it from global list
if not self._id is None:
MATERIAL_IDS.remove(self._id)
if material_id is None:
global AUTO_MATERIAL_ID
self._id = AUTO_MATERIAL_ID
MATERIAL_IDS.append(AUTO_MATERIAL_ID)
AUTO_MATERIAL_ID += 1
@ -243,7 +241,7 @@ class Material(object):
def get_all_nuclides(self):
nuclides = dict()
nuclides = {}
for nuclide_name, nuclide_tuple in self._nuclides.items():
nuclide = nuclide_tuple[0]
@ -262,13 +260,13 @@ class Material(object):
string += '{0: <16}{1}{2}'.format('\tDensity', '=\t', self._density)
string += ' [{0}]\n'.format(self._density_units)
string += '{0: <16}'.format('\tS(a,b) Tables') + '\n'
string += '{0: <16}\n'.format('\tS(a,b) Tables')
for sab in self._sab:
string += '{0: <16}{1}[{2}{3}]\n'.format('\tS(a,b)', '=\t',
sab[0], sab[1])
string += '{0: <16}'.format('\tNuclides') + '\n'
string += '{0: <16}\n'.format('\tNuclides')
for nuclide in self._nuclides:
percent = self._nuclides[nuclide][1]
@ -320,7 +318,7 @@ class Material(object):
def get_nuclides_xml(self, nuclides, distrib=False):
xml_elements = list()
xml_elements = []
for nuclide in nuclides.values():
xml_elements.append(self.get_nuclide_xml(nuclide, distrib))
@ -330,7 +328,7 @@ class Material(object):
def get_elements_xml(self, elements, distrib=False):
xml_elements = list()
xml_elements = []
for element in elements.values():
xml_elements.append(self.get_element_xml(element, distrib))
@ -365,7 +363,7 @@ class Material(object):
else:
subelement = ET.SubElement(element, "compositions")
comps = []
allnucs = self._nuclides.values() + self._elements.values()
dist_per_type = allnucs[0][2]
@ -416,7 +414,7 @@ class MaterialsFile(object):
def __init__(self):
# Initialize MaterialsFile class attributes
self._materials = list()
self._materials = []
self._default_xs = None
self._materials_file = ET.Element("materials")

9
src/utils/openmc/nuclide.py
View file

@ -1,5 +1,3 @@
#!/usr/bin/env python
from openmc.checkvalue import *
@ -38,10 +36,7 @@ class Nuclide(object):
def __hash__(self):
hashable = list()
hashable.append(self._name)
hashable.append(self._xs)
return hash(tuple(hashable))
return hash((self._name, self._xs))
def set_name(self, name):
@ -80,4 +75,4 @@ class Nuclide(object):
string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self._xs)
if self._zaid is not None:
string += '{0: <16}{1}{2}\n'.format('\tZAID', '=\t', self._zaid)
return string
return string

30
src/utils/openmc/opencg_compatible.py
View file

@ -1,60 +1,60 @@
#!/usr/bin/env python
import copy
import numpy as np
import opencg
import openmc
import opencg
import copy
import numpy as np
# A dictionary of all OpenMC Materials created
# Keys - Material IDs
# Values - Materials
OPENMC_MATERIALS = dict()
OPENMC_MATERIALS = {}
# A dictionary of all OpenCG Materials created
# Keys - Material IDs
# Values - Materials
OPENCG_MATERIALS = dict()
OPENCG_MATERIALS = {}
# A dictionary of all OpenMC Surfaces created
# Keys - Surface IDs
# Values - Surfaces
OPENMC_SURFACES = dict()
OPENMC_SURFACES = {}
# A dictionary of all OpenCG Surfaces created
# Keys - Surface IDs
# Values - Surfaces
OPENCG_SURFACES = dict()
OPENCG_SURFACES = {}
# A dictionary of all OpenMC Cells created
# Keys - Cell IDs
# Values - Cells
OPENMC_CELLS = dict()
OPENMC_CELLS = {}
# A dictionary of all OpenCG Cells created
# Keys - Cell IDs
# Values - Cells
OPENCG_CELLS = dict()
OPENCG_CELLS = {}
# A dictionary of all OpenMC Universes created
# Keys - Universes IDs
# Values - Universes
OPENMC_UNIVERSES = dict()
OPENMC_UNIVERSES = {}
# A dictionary of all OpenCG Universes created
# Keys - Universes IDs
# Values - Universes
OPENCG_UNIVERSES = dict()
OPENCG_UNIVERSES = {}
# A dictionary of all OpenMC Lattices created
# Keys - Lattice IDs
# Values - Lattices
OPENMC_LATTICES = dict()
OPENMC_LATTICES = {}
# A dictionary of all OpenCG Lattices created
# Keys - Lattice IDs
# Values - Lattices
OPENCG_LATTICES = dict()
OPENCG_LATTICES = {}
@ -410,7 +410,7 @@ def get_compatible_opencg_cells(opencg_cell, opencg_surface, halfspace):
raise ValueError(msg)
# Initialize an empty list for the new compatible cells
compatible_cells = list()
compatible_cells = []
# SquarePrism Surfaces
if opencg_surface._type in ['x-squareprism',

8
src/utils/openmc/plots.py
View file

@ -1,9 +1,9 @@
#!/usr/bin/env python
from xml.etree import ElementTree as ET
import numpy as np
from openmc.checkvalue import *
from openmc.clean_xml import *
from xml.etree import ElementTree as ET
import numpy as np
# A static variable for auto-generated Plot IDs
@ -400,7 +400,7 @@ class PlotsFile(object):
def __init__(self):
# Initialize PlotsFile class attributes
self._plots = list()
self._plots = []
self._plots_file = ET.Element("plots")

65
src/utils/openmc/settings.py
View file

@ -1,12 +1,11 @@
#!/usr/bin/env python
import collections
import warnings
from xml.etree import ElementTree as ET
import numpy as np
from openmc.checkvalue import *
from openmc.clean_xml import *
from xml.etree import ElementTree as ET
import numpy as np
import multiprocessing
class SettingsFile(object):
@ -154,16 +153,16 @@ class SettingsFile(object):
self._source_file = source_file
def set_source_space(self, type, params):
def set_source_space(self, stype, params):
if not is_string(type):
if not is_string(stype):
msg = 'Unable to set source space type to a non-string ' \
'value {0}'.format(type)
'value {0}'.format(stype)
raise ValueError(msg)
elif not type in ['box', 'point']:
elif not stype in ['box', 'point']:
msg = 'Unable to set source space type to {0} since it is not ' \
'box or point'.format(type)
'box or point'.format(stype)
raise ValueError(msg)
elif not isinstance(params, (tuple, list, np.ndarray)):
@ -183,20 +182,20 @@ class SettingsFile(object):
'is not an integer or floating point value'.format(param)
raise ValueError(msg)
self._source_space_type = type
self._source_space_type = stype
self._source_space_params = params
def set_source_angle(self, type, params=[]):
def set_source_angle(self, stype, params=[]):
if not is_string(type):
if not is_string(stype):
msg = 'Unable to set source angle type to a non-string ' \
'value {0}'.format(type)
'value {0}'.format(stype)
raise ValueError(msg)
elif not type in ['isotropic', 'monodirectional']:
elif not stype in ['isotropic', 'monodirectional']:
msg = 'Unable to set source angle type to {0} since it is not ' \
'isotropic or monodirectional'.format(type)
'isotropic or monodirectional'.format(stype)
raise ValueError(msg)
elif not isinstance(params, (tuple, list, np.ndarray)):
@ -204,12 +203,12 @@ class SettingsFile(object):
'not a Python list/tuple or NumPy array'.format(params)
raise ValueError(msg)
elif type is 'isotropic' and not params is None:
elif stype == 'isotropic' and not params is None:
msg = 'Unable to set source angle parameters since they are not ' \
'it is not supported for isotropic type sources'
raise ValueError(msg)
elif type is 'monodirectional' and len(params) != 3:
elif stype == 'monodirectional' and len(params) != 3:
msg = 'Unable to set source angle parameters to {0} ' \
'since 3 parameters are required for monodirectional ' \
'sources'.format(params)
@ -222,20 +221,20 @@ class SettingsFile(object):
'is not an integer or floating point value'.format(param)
raise ValueError(msg)
self._source_angle_type = type
self._source_angle_type = stype
self._source_angle_params = params
def set_source_energy(self, type, params=[]):
def set_source_energy(self, stype, params=[]):
if not is_string(type):
if not is_string(stype):
msg = 'Unable to set source energy type to a non-string ' \
'value {0}'.format(type)
'value {0}'.format(stype)
raise ValueError(msg)
elif not type in ['monoenergetic', 'watt', 'maxwell']:
elif not stype in ['monoenergetic', 'watt', 'maxwell']:
msg = 'Unable to set source energy type to {0} since it is not ' \
'monoenergetic, watt or maxwell'.format(type)
'monoenergetic, watt or maxwell'.format(stype)
raise ValueError(msg)
elif not isinstance(params, (tuple, list, np.ndarray)):
@ -243,19 +242,19 @@ class SettingsFile(object):
'is not a Python list/tuple or NumPy array'.format(params)
raise ValueError(msg)
elif type is 'monoenergetic' and not len(params) != 1:
elif stype == 'monoenergetic' and not len(params) != 1:
msg = 'Unable to set source energy parameters to {0} ' \
'since 1 paramater is required for monenergetic ' \
'sources'.format(params)
raise ValueError(msg)
elif type is 'watt' and len(params) != 2:
elif stype == 'watt' and len(params) != 2:
msg = 'Unable to set source energy parameters to {0} ' \
'since 2 parameters are required for monoenergetic ' \
'sources'.format(params)
raise ValueError(msg)
elif type is 'maxwell' and len(params) != 2:
elif stype == 'maxwell' and len(params) != 2:
msg = 'Unable to set source energy parameters to {0} since 1 ' \
'parameter is required for maxwell sources'.format(params)
raise ValueError(msg)
@ -268,7 +267,7 @@ class SettingsFile(object):
'value'.format(param)
raise ValueError(msg)
self._source_energy_type = type
self._source_energy_type = stype
self._source_energy_params = params
@ -437,9 +436,9 @@ class SettingsFile(object):
def set_energy_grid(self, energy_grid):
if not energy_grid in ['union', 'nuclide']:
if not energy_grid in ['nuclide', 'logarithm', 'material-union']:
msg = 'Unable to set energy grid to {0} which is neither ' \
'union nor nuclide'.format(energy_grid)
'nuclide, logarithm, nor material-union'.format(energy_grid)
raise ValueError(msg)
self._energy_grid = energy_grid
@ -808,7 +807,7 @@ class SettingsFile(object):
warnings.warn('This feature is not yet implemented in a release ' \
'version of openmc')
if not type(allow) == bool:
if not isinstance(allow, bool):
msg = 'Unable to set DD allow_leakage {0} which is ' \
'not a Python bool'.format(dimension)
raise ValueError(msg)
@ -822,7 +821,7 @@ class SettingsFile(object):
warnings.warn('This feature is not yet implemented in a release ' \
'version of openmc')
if not type(interactions) == bool:
if not isinstance(interactions, bool):
msg = 'Unable to set DD count_interactions {0} which is ' \
'not a Python bool'.format(dimension)
raise ValueError(msg)
@ -1206,7 +1205,7 @@ class SettingsFile(object):
subelement.text = 'true'
else:
subelement.text = 'false'
subelement = ET.SubElement(element, "count_interactions")
if self._dd_count_interactions:
subelement.text = 'true'

30
src/utils/openmc/statepoint.py
View file

@ -1,13 +1,13 @@
#!/usr/bin/env python
import copy
import struct
import struct, copy
import numpy as np
import scipy.stats
import openmc
from openmc.constants import *
class SourceSite(object):
def __init__(self):
@ -188,7 +188,7 @@ class StatePoint(object):
# Initialize dictionaries for the Meshes
# Keys - Mesh IDs
# Values - Mesh objects
self._meshes = dict()
self._meshes = {}
# Read the number of Meshes
self._n_meshes = self._get_int(path='tallies/meshes/n_meshes')[0]
@ -205,8 +205,8 @@ class StatePoint(object):
path='tallies/meshes/keys')
else:
self._mesh_keys = list()
self._mesh_ids = list()
self._mesh_keys = []
self._mesh_ids = []
# Build dictionary of Meshes
base = 'tallies/meshes/mesh '
@ -253,7 +253,7 @@ class StatePoint(object):
# Initialize dictionaries for the Tallies
# Keys - Tally IDs
# Values - Tally objects
self._tallies = dict()
self._tallies = {}
# Read the number of tallies
self._n_tallies = self._get_int(path='/tallies/n_tallies')[0]
@ -270,8 +270,8 @@ class StatePoint(object):
self._n_tallies, path='tallies/keys')
else:
self._tally_keys = list()
self._tally_ids = list()
self._tally_keys = []
self._tally_ids = []
base = 'tallies/tally '
@ -378,12 +378,12 @@ class StatePoint(object):
path='{0}{1}/n_user_score_bins'.format(base, tally_key))[0]
# Read scattering moment order strings (e.g., P3, Y-1,2, etc.)
moments = list()
moments = []
subbase = '{0}{1}/moments/'.format(base, tally_key)
# Extract the moment order string for each score
for k in range(len(scores)):
moment = self._get_string(8,
moment = self._get_string(8,
path='{0}order{1}'.format(subbase, k+1))
moment = moment.lstrip('[\'')
moment = moment.rstrip('\']')
@ -666,25 +666,25 @@ class StatePoint(object):
for filter in tally._filters:
if filter._type == 'surface':
surface_ids = list()
surface_ids = []
for bin in filter._bins:
surface_ids.append(summary.surfaces[bin]._id)
filter.set_bin_edges(surface_ids)
if filter._type in ['cell', 'distribcell']:
distribcell_ids = list()
distribcell_ids = []
for bin in filter._bins:
distribcell_ids.append(summary.cells[bin]._id)
filter.set_bin_edges(distribcell_ids)
if filter._type == 'universe':
universe_ids = list()
universe_ids = []
for bin in filter._bins:
universe_ids.append(summary.universes[bin]._id)
filter.set_bin_edges(universe_ids)
if filter._type == 'material':
material_ids = list()
material_ids = []
for bin in filter._bins:
material_ids.append(summary.materials[bin]._id)
filter.set_bin_edges(material_ids)

25
src/utils/openmc/summary.py
View file

@ -1,8 +1,7 @@
#!/usr/bin/env python
import numpy as np
import openmc
from openmc.opencg_compatible import get_opencg_geometry
import numpy as np
try:
import h5py
@ -65,7 +64,7 @@ class Summary(object):
# Initialize dictionary for each Nuclide
# Keys - Nuclide ZAIDs
# Values - Nuclide objects
self.nuclides = dict()
self.nuclides = {}
for key in self._f['nuclides'].keys():
@ -97,7 +96,7 @@ class Summary(object):
# Initialize dictionary for each Material
# Keys - Material keys
# Values - Material objects
self.materials = dict()
self.materials = {}
for key in self._f['materials'].keys():
@ -111,8 +110,8 @@ class Summary(object):
nuclides = self._f['materials'][key]['nuclides'][...]
n_sab = self._f['materials'][key]['n_sab'][0]
sab_names = list()
sab_xs = list()
sab_names = []
sab_xs = []
# Read the names of the S(a,b) tables for this Material
for i in range(1, n_sab+1):
@ -156,7 +155,7 @@ class Summary(object):
# Initialize dictionary for each Surface
# Keys - Surface keys
# Values - Surfacee objects
self.surfaces = dict()
self.surfaces = {}
for key in self._f['geometry/surfaces'].keys():
@ -239,7 +238,7 @@ class Summary(object):
# Initialize dictionary for each Cell
# Keys - Cell keys
# Values - Cell objects
self.cells = dict()
self.cells = {}
# Initialize dictionary for each Cell's fill
# (e.g., Material, Universe or Lattice ID)
@ -247,7 +246,7 @@ class Summary(object):
# the corresponding objects
# Keys - Cell keys
# Values - Filling Material, Universe or Lattice ID
self._cell_fills = dict()
self._cell_fills = {}
for key in self._f['geometry/cells'].keys():
@ -268,7 +267,7 @@ class Summary(object):
if 'surfaces' in self._f['geometry/cells'][key].keys():
surfaces = self._f['geometry/cells'][key]['surfaces'][...]
else:
surfaces = list()
surfaces = []
# Create this Cell
cell = openmc.Cell(cell_id=cell_id)
@ -310,7 +309,7 @@ class Summary(object):
# Initialize dictionary for each Universe
# Keys - Universe keys
# Values - Universe objects
self.universes = dict()
self.universes = {}
for key in self._f['geometry/universes'].keys():
@ -340,7 +339,7 @@ class Summary(object):
# Initialize lattices for each Lattice
# Keys - Lattice keys
# Values - Lattice objects
self.lattices = dict()
self.lattices = {}
for key in self._f['geometry/lattices'].keys():
@ -354,7 +353,7 @@ class Summary(object):
if lattice_type == 'rectangular':
dimension = self._f['geometry/lattices'][key]['n_cells'][...]
lower_left = \
self._f['geometry/lattices'][key]['lower_left'][...]
self._f['geometry/lattices'][key]['lower_left'][...]
pitch = self._f['geometry/lattices'][key]['pitch'][...]
outer = self._f['geometry/lattices'][key]['outer'][0]

9
src/utils/openmc/surface.py
View file

@ -1,8 +1,7 @@
#!/usr/bin/env python
from xml.etree import ElementTree as ET
from openmc.checkvalue import *
from openmc.constants import BC_TYPES
from xml.etree import ElementTree as ET
# A static variable for auto-generated Surface IDs
@ -27,11 +26,11 @@ class Surface(object):
# A dictionary of the quadratic surface coefficients
# Key - coefficeint name
# Value - coefficient value
self._coeffs = dict()
self._coeffs = {}
# An ordered list of the coefficient names to export to XML in the
# proper order
self._coeff_keys = list()
self._coeff_keys = []
self.set_id(surface_id)
self.set_boundary_type(bc_type)
@ -592,4 +591,4 @@ class ZCone(Cone):
# Initialize ZCone class attributes
super(ZCone, self).__init__(surface_id, bc_type, x0, y0, z0, R2, name=name)
self._type = 'z-cone'
self._type = 'z-cone'

39
src/utils/openmc/tallies.py
View file

@ -1,12 +1,13 @@
#!/usr/bin/env python
import copy
import os
from xml.etree import ElementTree as ET
import numpy as np
from openmc import Nuclide
from openmc.clean_xml import *
from openmc.checkvalue import *
from openmc.constants import *
from xml.etree import ElementTree as ET
import numpy as np
import os, copy
# "Static" variables for auto-generated Tally and Mesh IDs
@ -64,7 +65,7 @@ class Filter(object):
def __hash__(self):
hashable = list()
hashable = []
hashable.append(self._type)
hashable.append(self._bins)
return hash(tuple(hashable))
@ -532,9 +533,9 @@ class Tally(object):
# Initialize Tally class attributes
self._id = None
self._label = None
self._filters = list()
self._nuclides = list()
self._scores = list()
self._filters = []
self._nuclides = []
self._scores = []
self._estimator = None
self._num_score_bins = 0
@ -567,15 +568,15 @@ class Tally(object):
clone._mean = copy.deepcopy(self._mean, memo)
clone._std_dev = copy.deepcopy(self._std_dev, memo)
clone._filters = list()
clone._filters = []
for filter in self._filters:
clone.add_filter(copy.deepcopy(filter, memo))
clone._nuclides = list()
clone._nuclides = []
for nuclide in self._nuclides:
clone.add_nuclide(copy.deepcopy(nuclide, memo))
clone._scores = list()
clone._scores = []
for score in self._scores:
clone.add_score(score)
@ -612,7 +613,7 @@ class Tally(object):
def __hash__(self):
hashable = list()
hashable = []
for filter in self._filters:
hashable.append((filter._type, tuple(filter._bins)))
@ -1105,7 +1106,7 @@ class Tally(object):
tally_group.create_dataset('scores', data=np.array(self._scores))
# Add a string array of the nuclides to the HDF5 group
nuclides = list()
nuclides = []
for nuclide in self._nuclides:
nuclides.append(nuclide._name)
@ -1140,10 +1141,10 @@ class Tally(object):
if os.path.exists(filename) and append:
tally_results = pickle.load(file(filename, 'rb'))
else:
tally_results = dict()
tally_results = {}
# Create a nested dictionary within the file for this particular Tally
tally_results['Tally-{0}'.format(self._id)] = dict()
tally_results['Tally-{0}'.format(self._id)] = {}
tally_group = tally_results['Tally-{0}'.format(self._id)]
# Add basic Tally data to the nested dictionary
@ -1153,7 +1154,7 @@ class Tally(object):
tally_group['scores'] = np.array(self._scores)
# Add a string array of the nuclides to the HDF5 group
nuclides = list()
nuclides = []
for nuclide in self._nuclides:
nuclides.append(nuclide._name)
@ -1161,7 +1162,7 @@ class Tally(object):
tally_group['nuclides']= np.array(nuclides)
# Create a nested dictionary for the Filters
tally_group['filters'] = dict()
tally_group['filters'] = {}
filter_group = tally_group['filters']
for filter in self._filters:
@ -1182,8 +1183,8 @@ class TalliesFile(object):
def __init__(self):
# Initialize TalliesFile class attributes
self._tallies = list()
self._meshes = list()
self._tallies = []
self._meshes = []
self._tallies_file = ET.Element("tallies")

34
src/utils/openmc/universe.py
View file

@ -1,11 +1,11 @@
#!/usr/bin/env python
import abc
from collections import OrderedDict
from xml.etree import ElementTree as ET
import numpy as np
import openmc
from openmc.checkvalue import *
from xml.etree import ElementTree as ET
from collections import OrderedDict
import numpy as np
import abc
################################################################################
@ -33,7 +33,7 @@ class Cell(object):
self._name = None
self._fill = None
self._type = None
self._surfaces = dict()
self._surfaces = {}
self._rotation = None
self._translation = None
self._offset = None
@ -228,7 +228,7 @@ class Cell(object):
def get_all_nuclides(self):
nuclides = dict()
nuclides = {}
if self._type != 'void':
nuclides.update(self._fill.get_all_nuclides())
@ -238,7 +238,7 @@ class Cell(object):
def get_all_cells(self):
cells = dict()
cells = {}
if self._type == 'fill' or self._type == 'lattice':
cells.update(self._fill.get_all_cells())
@ -248,7 +248,7 @@ class Cell(object):
def get_all_universes(self):
universes = dict()
universes = {}
if self._type == 'fill':
universes[self._fill._id] = self._fill
@ -382,7 +382,7 @@ class Universe(object):
# Keys - Cell IDs
# Values - Cells
self._cells = dict()
self._cells = {}
# Keys - Cell IDs
# Values - Offsets
@ -485,7 +485,7 @@ class Universe(object):
def get_all_nuclides(self):
nuclides = dict()
nuclides = {}
# Append all Nuclides in each Cell in the Universe to the dictionary
for cell_id, cell in self._cells.items():
@ -496,7 +496,7 @@ class Universe(object):
def get_all_cells(self):
cells = dict()
cells = {}
# Add this Universe's cells to the dictionary
cells.update(self._cells)
@ -513,7 +513,7 @@ class Universe(object):
# Get all Cells in this Universe
cells = self.get_all_cells()
universes = dict()
universes = {}
# Append all Universes containing each Cell to the dictionary
for cell_id, cell in cells.items():
@ -635,7 +635,7 @@ class Lattice(object):
def get_unique_universes(self):
unique_universes = np.unique(self._universes.ravel())
universes = dict()
universes = {}
for universe in unique_universes:
universes[universe._id] = universe
@ -645,7 +645,7 @@ class Lattice(object):
def get_all_nuclides(self):
nuclides = dict()
nuclides = {}
# Get all unique Universes contained in each of the lattice cells
unique_universes = self.get_unique_universes()
@ -659,7 +659,7 @@ class Lattice(object):
def get_all_cells(self):
cells = dict()
cells = {}
unique_universes = self.get_unique_universes()
for universe_id, universe in unique_universes.items():
@ -672,7 +672,7 @@ class Lattice(object):
# Initialize a dictionary of all Universes contained by the Lattice
# in each nested Universe level
all_universes = dict()
all_universes = {}
# Get all unique Universes contained in each of the lattice cells
unique_universes = self.get_unique_universes()

4
tests/run_tests.py
View file

@ -179,6 +179,10 @@ class Test(object):
# Runs cmake when in non-script mode
def run_cmake(self):
os.environ['FC'] = self.fc
if self.petsc:
os.environ['PETSC_DIR'] = PETSC_DIR
if self.mpi:
os.environ['MPI_DIR'] = MPI_DIR
build_opts = self.build_opts.split()
self.cmake += build_opts
rc = call(self.cmake)

200
tests/test_score_total_yn/results_true.dat
View file

@ -60,6 +60,106 @@ tally 2:
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
8.795501E-01
1.600341E-01
-1.530661E-02
5.215322E-04
1.571095E-02
9.166057E-04
7.226989E-03
3.683499E-04
-2.456165E-02
1.571441E-04
1.387449E-02
5.600648E-04
-2.186870E-02
2.081164E-04
-1.106814E-02
1.065144E-04
-4.579593E-03
2.008225E-04
7.573253E-03
2.493075E-04
-1.505016E-02
1.451400E-04
1.503792E-02
9.539029E-05
-1.183668E-02
1.741393E-04
4.060912E-03
5.890591E-05
1.789747E-02
8.239749E-05
-2.168438E-02
1.555757E-04
-1.045826E-02
8.108402E-05
1.227413E-02
2.502871E-04
-2.806122E-02
1.886120E-04
-4.918718E-03
2.129038E-04
-1.014729E-02
6.914145E-05
-5.873760E-03
2.229738E-04
6.666510E-03
1.394147E-04
-3.245528E-03
1.550876E-04
-1.037659E-02
2.163837E-04
1.517577E+01
4.747271E+01
-2.463504E-01
@ -110,6 +210,56 @@ tally 2:
2.383142E-02
-1.463402E-01
1.526988E-02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.151504E+00
2.051857E+00
-4.923938E-02
@ -160,6 +310,56 @@ tally 2:
1.674419E-04
-1.682320E-02
8.389254E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
4.536316E+01
4.258781E+02
-6.747275E-01

1
tests/test_score_total_yn/tallies.xml
View file

@ -9,6 +9,7 @@
<tally id="2">
<filter type="cell" bins="10 21 22 23" />
<scores>total-y4</scores>
<nuclides>U-235 total</nuclides>
</tally>
</tallies>

8
tests/test_union_energy_grids/geometry.xml Normal file
View file

@ -0,0 +1,8 @@
<?xml version="1.0"?>
<geometry>
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" surfaces="-1" />
</geometry>

11
tests/test_union_energy_grids/materials.xml Normal file
View file

@ -0,0 +1,11 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U-235" xs="71c" ao="1.0" />
<nuclide name="H-1" xs="71c" ao="0.5" />
<nuclide name="C-Nat" xs="71c" ao="0.5" />
</material>
</materials>

25
tests/test_union_energy_grids/results.py Normal file
View file

@ -0,0 +1,25 @@
#!/usr/bin/env python
import sys
sys.path.insert(0, '../../src/utils')
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
print(sys.argv)
sp = StatePoint(sys.argv[1])
else:
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

2
tests/test_union_energy_grids/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
3.215828E-01 2.966835E-03

18
tests/test_union_energy_grids/settings.xml Normal file
View file

@ -0,0 +1,18 @@
<?xml version="1.0"?>
<settings>
<energy_grid>union</energy_grid>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
</settings>

59
tests/test_union_energy_grids/test_union_energy_grids.py Normal file
View file

@ -0,0 +1,59 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
import glob
from optparse import OptionParser
parser = OptionParser()
parser.add_option('--mpi_exec', dest='mpi_exec', default='')
parser.add_option('--mpi_np', dest='mpi_np', default='3')
parser.add_option('--exe', dest='exe')
(opts, args) = parser.parse_args()
cwd = os.getcwd()
def test_run():
if opts.mpi_exec != '':
proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE)
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0, 'OpenMC did not exit successfully.'
def test_created_statepoint():
statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*'))
assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.'
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\
'Statepoint file is not a binary or hdf5 file.'
def test_results():
statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*'))
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare, 'Results do not agree.'
def teardown():
output = glob.glob(os.path.join(cwd, 'statepoint.10.*'))
output.append(os.path.join(cwd, 'results_test.dat'))
for f in output:
if os.path.exists(f):
os.remove(f)
if __name__ == '__main__':
# test for openmc executable
if opts.exe is None:
raise Exception('Must specify OpenMC executable from command line with --exe.')
# run tests
try:
test_run()
test_created_statepoint()
test_results()
finally:
teardown()