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

Merge pull request #379 from wbinventor/squash-distribcells

Browse source 
Distributed cell (distribcell) tally filters
This commit is contained in:
Paul Romano 2015-05-29 08:56:51 +07:00
commit 3b189813aa
39 changed files with 1865 additions and 87 deletions
Download patch file Download diff file Expand all files Collapse all files

18
.gitignore vendored
View file

@ -39,11 +39,27 @@ src/templates/*.f90
# Test results error file
results_error.dat
# Test build files
tests/build/
tests/ctestscript.run
# HDF5 files
*.h5
# Build files
src/CMakeCache.txt
src/CMakeFiles/
src/bin/
src/cmake_install.cmake
src/install_manifest.txt
# Data downloaded from NNDC
data/nndc
#Images
*.ppm
# PyCharm project configuration files
.idea
.idea/*
.idea/*

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

@ -1218,8 +1218,9 @@ The ``<tally>`` element accepts the following sub-elements:
The ``filter`` element has the following attributes/sub-elements:
:type:
The type of the filter. Accepted options are "cell", "cellborn",
"material", "universe", "energy", "energyout", and "mesh".
The type of the filter. Accepted options are "cell", "cellborn",
"material", "universe", "energy", "energyout", "mesh", and
"distribcell".
:bins:
For each filter type, the corresponding ``bins`` entry is given as
@ -1258,6 +1259,15 @@ The ``<tally>`` element accepts the following sub-elements:
:mesh:
The ``id`` of a structured mesh to be tallied over.
:distribcell:
The single cell which should be tallied uniquely for all instances.
.. note::
The distribcell filter will take a single cell ID and will tally
each unique occurrence of that cell separately. This filter will
not accept more than one cell ID. It is not recommended to combine
this filter with a cell or mesh filter.
:nuclides:
If specified, the scores listed will be for particular nuclides, not the
summation of reactions from all nuclides. The format for nuclides should be

15
examples/python/lattice/hexagonal/build-xml.py
View file

@ -159,3 +159,18 @@ plot_file = openmc.PlotsFile()
plot_file.add_plot(plot_xy)
plot_file.add_plot(plot_yz)
plot_file.export_to_xml()
###############################################################################
# Exporting to OpenMC tallies.xml File
###############################################################################
# Instantiate a distribcell Tally
tally = openmc.Tally(tally_id=1)
tally.add_filter(openmc.Filter(type='distribcell', bins=[cell2.id]))
tally.add_score('total')
# Instantiate a TalliesFile, register Tally/Mesh, and export to XML
tallies_file = openmc.TalliesFile()
tallies_file.add_tally(tally)
tallies_file.export_to_xml()

19
src/constants.F90
View file

@ -297,16 +297,17 @@ module constants
integer, parameter :: NO_BIN_FOUND = -1
! Tally filter and map types
integer, parameter :: N_FILTER_TYPES = 8
integer, parameter :: N_FILTER_TYPES = 9
integer, parameter :: &
FILTER_UNIVERSE = 1, &
FILTER_MATERIAL = 2, &
FILTER_CELL = 3, &
FILTER_CELLBORN = 4, &
FILTER_SURFACE = 5, &
FILTER_MESH = 6, &
FILTER_ENERGYIN = 7, &
FILTER_ENERGYOUT = 8
FILTER_UNIVERSE = 1, &
FILTER_MATERIAL = 2, &
FILTER_CELL = 3, &
FILTER_CELLBORN = 4, &
FILTER_SURFACE = 5, &
FILTER_MESH = 6, &
FILTER_ENERGYIN = 7, &
FILTER_ENERGYOUT = 8, &
FILTER_DISTRIBCELL = 9
! Tally surface current directions
integer, parameter :: &

347
src/geometry.F90
View file

@ -198,6 +198,12 @@ contains
p % coord => p % coord % next
p % coord % universe = c % fill
! Determine all distribcell offsets for this cell level
if (.not. associated(p % coord % mapping)) then
allocate(p % coord % mapping(n_maps))
end if
p % coord % mapping(:) = c % offset(:)
! Apply translation
if (allocated(c % translation)) then
p % coord % xyz = p % coord % xyz - c % translation
@ -253,6 +259,15 @@ contains
! Move particle to next level and search for the lower cells.
p % coord => p % coord % next
! Determine all distribcell offsets for this cell level
if (lat % are_valid_indices(i_xyz)) then
if (.not. associated(p % coord % mapping)) then
allocate(p % coord % mapping(n_maps))
end if
p % coord % mapping(:) = lat % offset(:, i_xyz(1), i_xyz(2), i_xyz(3))
end if
call find_cell(p, found)
if (.not. found) exit
@ -612,8 +627,16 @@ contains
end if
else OUTSIDE_LAT
! Find cell in next lattice element
p % coord % universe = lat % universes(i_xyz(1), i_xyz(2), i_xyz(3))
! Determine all distribcell offsets for this lattice cell
if (.not. associated(p % coord % mapping)) then
allocate(p % coord % mapping(n_maps))
end if
p % coord % mapping(:) = lat % offset(:, i_xyz(1), i_xyz(2), i_xyz(3))
call find_cell(p, found)
if (.not. found) then
! In some circumstances, a particle crossing the corner of a cell may
@ -1582,5 +1605,329 @@ contains
end if
end subroutine handle_lost_particle
!===============================================================================
! CALC_OFFSETS calculates and stores the offsets in all fill cells. This
! routine is called once upon initialization.
!===============================================================================
subroutine calc_offsets(goal, map, univ, counts, found)
integer, intent(inout) :: goal ! target universe ID
integer, intent(inout) :: map ! map index in vector of maps
type(Universe), intent(in) :: univ ! universe searching in
integer, intent(inout) :: counts(:,:) ! target count
logical, intent(inout) :: found(:,:) ! target found
integer :: i ! index over cells
integer :: j, k, m ! indices in lattice
integer :: n ! number of cells to search
integer :: offset ! total offset for a given cell
integer :: cell_index ! index in cells array
type(Cell), pointer :: c ! pointer to current cell
type(Universe), pointer :: next_univ ! next universe to cycle through
class(Lattice), pointer :: lat ! pointer to current lattice
n = univ % n_cells
offset = 0
do i = 1, n
cell_index = univ % cells(i)
! get pointer to cell
c => cells(cell_index)
! ====================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH
if (c % type == CELL_NORMAL) then
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
elseif (c % type == CELL_FILL) then
! Set offset for the cell on this level
c % offset(map) = offset
! Count contents of this cell
next_univ => universes(c % fill)
offset = offset + count_target(next_univ, counts, found, goal, map)
! Move into the next universe
next_univ => universes(c % fill)
c => cells(cell_index)
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == CELL_LATTICE) then
! Set current lattice
lat => lattices(c % fill) % obj
select type (lat)
type is (RectLattice)
! Loop over lattice coordinates
do j = 1, lat % n_cells(1)
do k = 1, lat % n_cells(2)
do m = 1, lat % n_cells(3)
lat % offset(map, j, k, m) = offset
next_univ => universes(lat % universes(j, k, m))
offset = offset + &
count_target(next_univ, counts, found, goal, map)
end do
end do
end do
type is (HexLattice)
! Loop over lattice coordinates
do m = 1, lat % n_axial
do k = 1, 2*lat % n_rings - 1
do j = 1, 2*lat % n_rings - 1
! This array location is never used
if (j + k < lat % n_rings + 1) then
cycle
! This array location is never used
else if (j + k > 3*lat % n_rings - 1) then
cycle
else
lat % offset(map, j, k, m) = offset
next_univ => universes(lat % universes(j, k, m))
offset = offset + &
count_target(next_univ, counts, found, goal, map)
end if
end do
end do
end do
end select
end if
end do
end subroutine calc_offsets
!===============================================================================
! COUNT_TARGET recursively totals the numbers of occurances of a given
! universe ID beginning with the universe given.
!===============================================================================
recursive function count_target(univ, counts, found, goal, map) result(count)
type(Universe), intent(inout) :: univ ! universe to search through
integer, intent(inout) :: counts(:,:) ! target count
logical, intent(inout) :: found(:,:) ! target found
integer, intent(inout) :: goal ! target universe ID
integer, intent(inout) :: map ! current map
integer :: i ! index over cells
integer :: j, k, m ! indices in lattice
integer :: n ! number of cells to search
integer :: cell_index ! index in cells array
integer :: count ! number of times target located
type(Cell), pointer :: c ! pointer to current cell
type(Universe), pointer :: next_univ ! next univ to loop through
class(Lattice), pointer :: lat ! pointer to current lattice
! Don't research places already checked
if (found(universe_dict % get_key(univ % id), map)) then
count = counts(universe_dict % get_key(univ % id), map)
return
end if
! If this is the target, it can't contain itself.
! Count = 1, then quit
if (univ % id == goal) then
count = 1
counts(universe_dict % get_key(univ % id), map) = 1
found(universe_dict % get_key(univ % id), map) = .true.
return
end if
count = 0
n = univ % n_cells
do i = 1, n
cell_index = univ % cells(i)
! get pointer to cell
c => cells(cell_index)
! ====================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH
if (c % type == CELL_NORMAL) then
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
elseif (c % type == CELL_FILL) then
next_univ => universes(c % fill)
! Found target - stop since target cannot contain itself
if (next_univ % id == goal) then
count = count + 1
return
end if
count = count + count_target(next_univ, counts, found, goal, map)
c => cells(cell_index)
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == CELL_LATTICE) then
! Set current lattice
lat => lattices(c % fill) % obj
select type (lat)
type is (RectLattice)
! Loop over lattice coordinates
do j = 1, lat % n_cells(1)
do k = 1, lat % n_cells(2)
do m = 1, lat % n_cells(3)
next_univ => universes(lat % universes(j, k, m))
! Found target - stop since target cannot contain itself
if (next_univ % id == goal) then
count = count + 1
cycle
end if
count = count + &
count_target(next_univ, counts, found, goal, map)
end do
end do
end do
type is (HexLattice)
! Loop over lattice coordinates
do m = 1, lat % n_axial
do k = 1, 2*lat % n_rings - 1
do j = 1, 2*lat % n_rings - 1
! This array location is never used
if (j + k < lat % n_rings + 1) then
cycle
! This array location is never used
else if (j + k > 3*lat % n_rings - 1) then
cycle
else
next_univ => universes(lat % universes(j, k, m))
! Found target - stop since target cannot contain itself
if (next_univ % id == goal) then
count = count + 1
cycle
end if
count = count + &
count_target(next_univ, counts, found, goal, map)
end if
end do
end do
end do
end select
end if
end do
counts(universe_dict % get_key(univ % id), map) = count
found(universe_dict % get_key(univ % id), map) = .true.
end function count_target
!===============================================================================
! COUNT_INSTANCE recursively totals the number of occurrences of all cells
! beginning with the universe given.
!===============================================================================
recursive subroutine count_instance(univ)
type(Universe), intent(in) :: univ ! universe to search through
integer :: i ! index over cells
integer :: j, k, m ! indices in lattice
integer :: n ! number of cells to search
integer :: cell_index ! index in cells array
type(Cell), pointer :: c ! pointer to current cell
type(Universe), pointer :: next_univ ! next universe to loop through
class(Lattice), pointer :: lat ! pointer to current lattice
n = univ % n_cells
do i = 1, n
cell_index = univ % cells(i)
! get pointer to cell
c => cells(cell_index)
c % instances = c % instances + 1
! ====================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH
if (c % type == CELL_NORMAL) then
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
elseif (c % type == CELL_FILL) then
next_univ => universes(c % fill)
call count_instance(next_univ)
c => cells(cell_index)
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == CELL_LATTICE) then
! Set current lattice
lat => lattices(c % fill) % obj
select type (lat)
type is (RectLattice)
! Loop over lattice coordinates
do j = 1, lat % n_cells(1)
do k = 1, lat % n_cells(2)
do m = 1, lat % n_cells(3)
next_univ => universes(lat % universes(j, k, m))
call count_instance(next_univ)
end do
end do
end do
type is (HexLattice)
! Loop over lattice coordinates
do m = 1, lat % n_axial
do k = 1, 2*lat % n_rings - 1
do j = 1, 2*lat % n_rings - 1
! This array location is never used
if (j + k < lat % n_rings + 1) then
cycle
! This array location is never used
else if (j + k > 3*lat % n_rings - 1) then
cycle
else
next_univ => universes(lat % universes(j, k, m))
call count_instance(next_univ)
end if
end do
end do
end do
end select
end if
end do
end subroutine count_instance
end module geometry

35
src/geometry_header.F90
View file

@ -7,13 +7,13 @@ module geometry_header
!===============================================================================
type Universe
integer :: id ! Unique ID
integer :: type ! Type
integer :: n_cells ! # of cells within
integer, allocatable :: cells(:) ! List of cells within
real(8) :: x0 ! Translation in x-coordinate
real(8) :: y0 ! Translation in y-coordinate
real(8) :: z0 ! Translation in z-coordinate
integer :: id ! Unique ID
integer :: type ! Type
integer :: n_cells ! # of cells within
integer, allocatable :: cells(:) ! List of cells within
real(8) :: x0 ! Translation in x-coordinate
real(8) :: y0 ! Translation in y-coordinate
real(8) :: z0 ! Translation in z-coordinate
end type Universe
!===============================================================================
@ -28,6 +28,7 @@ module geometry_header
integer :: outside ! Material to fill area outside
integer :: outer ! universe to tile outside the lat
logical :: is_3d ! Lattice has cells on z axis
integer, allocatable :: offset(:,:,:,:) ! Distribcell offsets
contains
@ -131,17 +132,19 @@ module geometry_header
!===============================================================================
type Cell
integer :: id ! Unique ID
integer :: id ! Unique ID
character(len=52) :: name = "" ! User-defined name
integer :: type ! Type of cell (normal, universe, lattice)
integer :: universe ! universe # this cell is in
integer :: fill ! universe # filling this cell
integer :: material ! Material within cell (0 for universe)
integer :: n_surfaces ! Number of surfaces within
integer :: type ! Type of cell (normal, universe, lattice)
integer :: universe ! universe # this cell is in
integer :: fill ! universe # filling this cell
integer :: instances ! number of instances of this cell in the geom
integer :: material ! Material within cell (0 for universe)
integer :: n_surfaces ! Number of surfaces within
integer, allocatable :: offset (:) ! Distribcell offset for tally counter
integer, allocatable :: &
& surfaces(:) ! List of surfaces bounding cell -- note that
! parentheses, union, etc operators will be listed
! here too
& surfaces(:) ! List of surfaces bounding cell -- note that
! parentheses, union, etc operators will be listed
! here too
! Rotation matrix and translation vector
real(8), allocatable :: translation(:)

3
src/global.F90
View file

@ -302,6 +302,9 @@ module global
! Particle restart run
logical :: particle_restart_run = .false.
! Number of distribcell maps
integer :: n_maps
! Write out initial source
logical :: write_initial_source = .false.

35
src/hdf5_summary.F90
View file

@ -154,6 +154,14 @@ contains
call su % write_data(universes(c % fill) % id, "fill", &
group="geometry/cells/cell " // trim(to_str(c % id)))
call su % write_data(size(c % offset), "maps", &
group="geometry/cells/cell " // trim(to_str(c % id)))
if (size(c % offset) > 0) then
call su % write_data(c % offset, "offset", &
length=size(c % offset), &
group="geometry/cells/cell " // trim(to_str(c % id)))
end if
if (allocated(c % translation)) then
call su % write_data(1, "translated", &
group="geometry/cells/cell " // trim(to_str(c % id)))
@ -164,7 +172,7 @@ contains
group="geometry/cells/cell " // trim(to_str(c % id)))
end if
if (allocated(c % rotation_matrix)) then
if (allocated(c % rotation)) then
call su % write_data(1, "rotated", &
group="geometry/cells/cell " // trim(to_str(c % id)))
call su % write_data(c % rotation, "rotation", length=3, &
@ -352,6 +360,16 @@ contains
call su % write_data(lat % outer, "outer", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
call su % write_data(size(lat % offset), "offset_size", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
call su % write_data(size(lat % offset,1), "maps", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
if (size(lat % offset) > 0) then
call su % write_data(lat % offset, "offsets", &
length=shape(lat % offset), &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
end if
! Write lattice universes.
allocate(lattice_universes(lat % n_cells(1), lat % n_cells(2), &
@ -398,13 +416,23 @@ contains
call su % write_data(lat % outer, "outer", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
call su % write_data(size(lat % offset), "offset_size", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
call su % write_data(size(lat % offset,1), "maps", &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
if (size(lat % offset) > 0) then
call su % write_data(lat % offset, "offsets", &
length=shape(lat % offset), &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
end if
! Write lattice universes.
allocate(lattice_universes(2*lat % n_rings - 1, 2*lat % n_rings - 1, &
&lat % n_axial))
do j = 1, lat % n_axial
do m = 1, lat % n_axial
do k = 1, 2*lat % n_rings - 1
do m = 1, 2*lat % n_rings - 1
do j = 1, 2*lat % n_rings - 1
if (j + k < lat % n_rings + 1) then
! This array position is never used; put a -1 to indicate this
lattice_universes(j,k,m) = -1
@ -488,6 +516,7 @@ contains
call su % write_data(m % i_sab_tables, "i_sab_tables", &
length=m % n_sab, &
group="materials/material " // trim(to_str(m % id)))
do j = 1, m % n_sab
call su % write_data(m % sab_names(j), to_str(j), &
group="materials/material " // &

227
src/initialize.F90
View file

@ -4,9 +4,10 @@ module initialize
use bank_header, only: Bank
use constants
use dict_header, only: DictIntInt, ElemKeyValueII
use set_header, only: SetInt
use energy_grid, only: logarithmic_grid, grid_method, unionized_grid
use error, only: fatal_error, warning
use geometry, only: neighbor_lists
use geometry, only: neighbor_lists, count_instance, calc_offsets
use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,&
&BASE_UNIVERSE
use global
@ -21,7 +22,7 @@ module initialize
use source, only: initialize_source
use state_point, only: load_state_point
use string, only: to_str, str_to_int, starts_with, ends_with
use tally_header, only: TallyObject, TallyResult
use tally_header, only: TallyObject, TallyResult, TallyFilter
use tally_initialize, only: configure_tallies
#ifdef MPI
@ -92,6 +93,9 @@ contains
! Use dictionaries to redefine index pointers
call adjust_indices()
! Initialize distribcell_filters
call prepare_distribcell()
! After reading input and basic geometry setup is complete, build lists of
! neighboring cells for efficient tracking
call neighbor_lists()
@ -688,6 +692,17 @@ contains
FILTER_LOOP: do j = 1, t % n_filters
select case (t % filters(j) % type)
case (FILTER_DISTRIBCELL)
do k = 1, size(t % filters(j) % int_bins)
id = t % filters(j) % int_bins(k)
if (cell_dict % has_key(id)) then
t % filters(j) % int_bins(k) = cell_dict % get_key(id)
else
call fatal_error("Could not find cell " // trim(to_str(id)) // &
" specified on tally " // trim(to_str(t % id)))
end if
end do
case (FILTER_CELL, FILTER_CELLBORN)
do k = 1, t % filters(j) % n_bins
@ -902,4 +917,212 @@ contains
end subroutine allocate_banks
!===============================================================================
! PREPARE_DISTRIBCELL initializes any distribcell filters present and sets the
! offsets for distribcells
!===============================================================================
subroutine prepare_distribcell()
integer :: i, j ! Tally, filter loop counters
integer :: n_filt ! Number of filters originally in tally
logical :: count_all ! Count all cells
type(TallyObject), pointer :: tally ! Current tally
type(Universe), pointer :: univ ! Pointer to universe
type(Cell), pointer :: c ! Pointer to cell
integer, allocatable :: univ_list(:) ! Target offsets
integer, allocatable :: counts(:,:) ! Target count
logical, allocatable :: found(:,:) ! Target found
count_all = .false.
! Loop over tallies
do i = 1, n_tallies
! Get pointer to tally
tally => tallies(i)
n_filt = tally % n_filters
! Loop over the filters to determine how many additional filters
! need to be added to this tally
do j = 1, tally % n_filters
! Determine type of filter
if (tally % filters(j) % type == FILTER_DISTRIBCELL) then
count_all = .true.
if (size(tally % filters(j) % int_bins) > 1) then
call fatal_error("A distribcell filter was specified with &
&multiple bins. This feature is not supported.")
end if
end if
end do
end do
if (count_all) then
univ => universes(BASE_UNIVERSE)
! sum the number of occurrences of all cells
call count_instance(univ)
! Loop over tallies
do i = 1, n_tallies
! Get pointer to tally
tally => tallies(i)
! Initialize the filters
do j = 1, tally % n_filters
! Set the number of bins to the number of instances of the cell
if (tally % filters(j) % type == FILTER_DISTRIBCELL) then
c => cells(tally % filters(j) % int_bins(1))
tally % filters(j) % n_bins = c % instances
end if
end do
end do
end if
! Allocate offset maps at each level in the geometry
call allocate_offsets(univ_list, counts, found)
! Calculate offsets for each target distribcell
do i = 1, n_maps
do j = 1, n_universes
univ => universes(j)
call calc_offsets(univ_list(i), i, univ, counts, found)
end do
end do
! Deallocate temporary target variable arrays
deallocate(counts)
deallocate(found)
deallocate(univ_list)
end subroutine prepare_distribcell
!===============================================================================
! ALLOCATE_OFFSETS determines the number of maps needed and allocates required
! memory for distribcell offset tables
!===============================================================================
recursive subroutine allocate_offsets(univ_list, counts, found)
integer, intent(out), allocatable :: univ_list(:) ! Target offsets
integer, intent(out), allocatable :: counts(:,:) ! Target count
logical, intent(out), allocatable :: found(:,:) ! Target found
integer :: i, j, k, l, m ! Loop counters
type(SetInt) :: cell_list ! distribells to track
type(Universe), pointer :: univ ! pointer to universe
class(Lattice), pointer :: lat ! pointer to lattice
type(TallyObject), pointer :: tally ! pointer to tally
type(TallyFilter), pointer :: filter ! pointer to filter
! Begin gathering list of cells in distribcell tallies
n_maps = 0
! Populate list of distribcells to track
do i = 1, n_tallies
tally => tallies(i)
do j = 1, tally % n_filters
filter => tally % filters(j)
if (filter % type == FILTER_DISTRIBCELL) then
if (.not. cell_list % contains(filter % int_bins(1))) then
call cell_list % add(filter % int_bins(1))
end if
end if
end do
end do
! Compute the number of unique universes containing these distribcells
! to determine the number of offset tables to allocate
do i = 1, n_universes
univ => universes(i)
do j = 1, univ % n_cells
if (cell_list % contains(univ % cells(j))) then
n_maps = n_maps + 1
end if
end do
end do
! Allocate the list of offset tables for each unique universe
allocate(univ_list(n_maps))
! Allocate list to accumulate target distribccell counts in each universe
allocate(counts(n_universes, n_maps))
! Allocate list to track if target distribcells are found in each universe
allocate(found(n_universes, n_maps))
counts(:,:) = 0
found(:,:) = .false.
k = 1
do i = 1, n_universes
univ => universes(i)
do j = 1, univ % n_cells
if (cell_list % contains(univ % cells(j))) then
! Loop over all tallies
do l = 1, n_tallies
tally => tallies(l)
do m = 1, tally % n_filters
filter => tally % filters(m)
! Loop over only distribcell filters
! If filter points to cell we just found, set offset index
if (filter % type == FILTER_DISTRIBCELL) then
if (filter % int_bins(1) == univ % cells(j)) then
filter % offset = k
end if
end if
end do
end do
univ_list(k) = univ % id
k = k + 1
end if
end do
end do
! Allocate the offset tables for lattices
do i = 1, n_lattices
lat => lattices(i) % obj
select type(lat)
type is (RectLattice)
allocate(lat % offset(n_maps, lat % n_cells(1), lat % n_cells(2), &
lat % n_cells(3)))
type is (HexLattice)
allocate(lat % offset(n_maps, 2 * lat % n_rings - 1, &
2 * lat % n_rings - 1, lat % n_axial))
end select
lat % offset(:, :, :, :) = 0
end do
! Allocate offset table for fill cells
do i = 1, n_cells
if (cells(i) % material == NONE) then
allocate(cells(i) % offset(n_maps))
end if
end do
end subroutine allocate_offsets
end module initialize

26
src/input_xml.F90
View file

@ -986,7 +986,6 @@ contains
integer :: n_cells_in_univ
integer :: coeffs_reqd
integer, allocatable :: temp_int_array(:)
real(8) :: temp_double_array3(3)
real(8) :: phi, theta, psi
logical :: file_exists
logical :: boundary_exists
@ -1044,6 +1043,9 @@ contains
do i = 1, n_cells
c => cells(i)
! Initialize the number of cell instances - this is a base case for distribcells
c % instances = 0
! Get pointer to i-th cell node
call get_list_item(node_cell_list, i, node_cell)
@ -1140,11 +1142,11 @@ contains
end if
! Copy rotation angles in x,y,z directions
call get_node_array(node_cell, "rotation", temp_double_array3)
c % rotation = temp_double_array3
phi = -temp_double_array3(1) * PI/180.0_8
theta = -temp_double_array3(2) * PI/180.0_8
psi = -temp_double_array3(3) * PI/180.0_8
allocate(c % rotation(3))
call get_node_array(node_cell, "rotation", c % rotation)
phi = -c % rotation(1) * PI/180.0_8
theta = -c % rotation(2) * PI/180.0_8
psi = -c % rotation(3) * PI/180.0_8
! Calculate rotation matrix based on angles given
allocate(c % rotation_matrix(3,3))
@ -2361,6 +2363,18 @@ contains
! Determine type of filter
select case (temp_str)
case ('distribcell')
! Set type of filter
t % filters(j) % type = FILTER_DISTRIBCELL
! Going to add new filters to this tally if n_words > 1
! Allocate and store bins
allocate(t % filters(j) % int_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % int_bins)
case ('cell')
! Set type of filter
t % filters(j) % type = FILTER_CELL

293
src/output.F90
View file

@ -5,7 +5,7 @@ module output
use ace_header, only: Nuclide, Reaction, UrrData
use constants
use endf, only: reaction_name
use error, only: warning
use error, only: fatal_error, warning
use geometry_header, only: Cell, Universe, Surface, Lattice, RectLattice, &
&HexLattice, BASE_UNIVERSE
use global
@ -754,6 +754,16 @@ contains
write(unit_,*) ' Estimator: Track-length'
end select
! Write any cells bins if present
j = t % find_filter(FILTER_DISTRIBCELL)
if (j > 0) then
string = ""
id = t % filters(j) % int_bins(1)
c => cells(id)
string = trim(string) // ' ' // trim(to_str(c % id))
write(unit_, *) ' Distribcell Bins:' // trim(string)
end if
! Write any cells bins if present
j = t % find_filter(FILTER_CELL)
if (j > 0) then
@ -863,7 +873,6 @@ contains
end do
write(unit_,*)
! Write score bins
string = ""
j = 0
@ -1712,7 +1721,7 @@ contains
real(8) :: t_value ! t-values for confidence intervals
real(8) :: alpha ! significance level for CI
character(MAX_FILE_LEN) :: filename ! name of output file
character(15) :: filter_name(N_FILTER_TYPES) ! names of tally filters
character(16) :: filter_name(N_FILTER_TYPES) ! names of tally filters
character(36) :: score_names(N_SCORE_TYPES) ! names of scoring function
character(36) :: score_name ! names of scoring function
! to be applied at write-time
@ -1722,14 +1731,15 @@ contains
if (n_tallies == 0) return
! Initialize names for tally filter types
filter_name(FILTER_UNIVERSE) = "Universe"
filter_name(FILTER_MATERIAL) = "Material"
filter_name(FILTER_CELL) = "Cell"
filter_name(FILTER_CELLBORN) = "Birth Cell"
filter_name(FILTER_SURFACE) = "Surface"
filter_name(FILTER_MESH) = "Mesh"
filter_name(FILTER_ENERGYIN) = "Incoming Energy"
filter_name(FILTER_ENERGYOUT) = "Outgoing Energy"
filter_name(FILTER_UNIVERSE) = "Universe"
filter_name(FILTER_MATERIAL) = "Material"
filter_name(FILTER_DISTRIBCELL) = "Distributed Cell"
filter_name(FILTER_CELL) = "Cell"
filter_name(FILTER_CELLBORN) = "Birth Cell"
filter_name(FILTER_SURFACE) = "Surface"
filter_name(FILTER_MESH) = "Mesh"
filter_name(FILTER_ENERGYIN) = "Incoming Energy"
filter_name(FILTER_ENERGYOUT) = "Outgoing Energy"
! Initialize names for scores
score_names(abs(SCORE_FLUX)) = "Flux"
@ -2120,14 +2130,16 @@ contains
type(TallyObject), pointer :: t ! tally object
integer, intent(in) :: i_filter ! index in filters array
character(30) :: label ! user-specified identifier
character(100) :: label ! user-specified identifier
integer :: i ! index in cells/surfaces/etc array
integer :: bin
integer :: offset
integer, allocatable :: ijk(:) ! indices in mesh
real(8) :: E0 ! lower bound for energy bin
real(8) :: E1 ! upper bound for energy bin
type(StructuredMesh), pointer :: m => null()
type(Universe), pointer :: univ => null()
bin = matching_bins(i_filter)
@ -2141,6 +2153,13 @@ contains
case (FILTER_CELL, FILTER_CELLBORN)
i = t % filters(i_filter) % int_bins(bin)
label = to_str(cells(i) % id)
case (FILTER_DISTRIBCELL)
label = ''
univ => universes(BASE_UNIVERSE)
offset = 0
call find_offset(t % filters(i_filter) % offset, &
t % filters(i_filter) % int_bins(1), &
univ, bin-1, offset, label)
case (FILTER_SURFACE)
i = t % filters(i_filter) % int_bins(bin)
label = to_str(surfaces(i) % id)
@ -2162,5 +2181,255 @@ contains
end select
end function get_label
!===============================================================================
! FIND_OFFSET uses a given map number, a target cell ID, and a target offset
! to build a string which is the path from the base universe to the target cell
! with the given offset
!===============================================================================
recursive subroutine find_offset(map, goal, univ, final, offset, path)
integer, intent(in) :: map ! Index in maps vector
integer, intent(in) :: goal ! The target cell ID
type(Universe), pointer, intent(in) :: univ ! Universe to begin search
integer, intent(in) :: final ! Target offset
integer, intent(inout) :: offset ! Current offset
character(100) :: path ! Path to offset
integer :: i, j ! Index over cells
integer :: k, l, m ! Indices in lattice
integer :: old_k, old_l, old_m ! Previous indices in lattice
integer :: n_x, n_y, n_z ! Lattice cell array dimensions
integer :: n ! Number of cells to search
integer :: cell_index ! Index in cells array
integer :: lat_offset ! Offset from lattice
integer :: temp_offset ! Looped sum of offsets
logical :: this_cell = .false. ! Advance in this cell?
logical :: later_cell = .false. ! Fill cells after this one?
type(Cell), pointer:: c ! Pointer to current cell
type(Universe), pointer :: next_univ ! Next universe to loop through
class(Lattice), pointer :: lat ! Pointer to current lattice
n = univ % n_cells
! Write to the geometry stack
if (univ%id == 0) then
path = trim(path) // to_str(univ%id)
else
path = trim(path) // "->" // to_str(univ%id)
end if
! Look through all cells in this universe
do i = 1, n
cell_index = univ % cells(i)
c => cells(cell_index)
! If the cell ID matches the goal and the offset matches final,
! write to the geometry stack
if (cell_dict % get_key(c % id) == goal .AND. offset == final) then
path = trim(path) // "->" // to_str(c%id)
return
end if
end do
! Find the fill cell or lattice cell that we need to enter
do i = 1, n
later_cell = .false.
cell_index = univ % cells(i)
c => cells(cell_index)
this_cell = .false.
! If we got here, we still think the target is in this universe
! or further down, but it's not this exact cell.
! Compare offset to next cell to see if we should enter this cell
if (i /= n) then
do j = i+1, n
cell_index = univ % cells(j)
c => cells(cell_index)
! Skip normal cells which do not have offsets
if (c % type == CELL_NORMAL) then
cycle
end if
! Break loop once we've found the next cell with an offset
exit
end do
! Ensure we didn't just end the loop by iteration
if (c % type /= CELL_NORMAL) then
! There are more cells in this universe that it could be in
later_cell = .true.
! Two cases, lattice or fill cell
if (c % type == CELL_FILL) then
temp_offset = c % offset(map)
! Get the offset of the first lattice location
else
lat => lattices(c % fill) % obj
temp_offset = lat % offset(map, 1, 1, 1)
end if
! If the final offset is in the range of offset - temp_offset+offset
! then the goal is in this cell
if (final < temp_offset + offset) then
this_cell = .true.
end if
end if
end if
if (n == 1 .and. c % type /= CELL_NORMAL) then
this_cell = .true.
end if
if (.not. later_cell) then
this_cell = .true.
end if
! Get pointer to THIS cell because target must be in this cell
if (this_cell) then
cell_index = univ % cells(i)
c => cells(cell_index)
path = trim(path) // "->" // to_str(c%id)
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
if (c % type == CELL_FILL) then
! Enter this cell to update the current offset
offset = c % offset(map) + offset
next_univ => universes(c % fill)
call find_offset(map, goal, next_univ, final, offset, path)
return
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type == CELL_LATTICE) then
! Set current lattice
lat => lattices(c % fill) % obj
select type (lat)
! ==================================================================
! RECTANGULAR LATTICES
type is (RectLattice)
! Write to the geometry stack
path = trim(path) // "->" // to_str(lat%id)
n_x = lat % n_cells(1)
n_y = lat % n_cells(2)
n_z = lat % n_cells(3)
old_m = 1
old_l = 1
old_k = 1
! Loop over lattice coordinates
do k = 1, n_x
do l = 1, n_y
do m = 1, n_z
if (final >= lat % offset(map, k, l, m) + offset) then
if (k == n_x .and. l == n_y .and. m == n_z) then
! This is last lattice cell, so target must be here
lat_offset = lat % offset(map, k, l, m)
offset = offset + lat_offset
next_univ => universes(lat % universes(k, l, m))
path = trim(path) // "(" // trim(to_str(k)) // &
"," // trim(to_str(l)) // "," // &
trim(to_str(m)) // ")"
call find_offset(map, goal, next_univ, final, offset, path)
return
else
old_m = m
old_l = l
old_k = k
cycle
end if
else
! Target is at this lattice position
lat_offset = lat % offset(map, old_k, old_l, old_m)
offset = offset + lat_offset
next_univ => universes(lat % universes(old_k, old_l, old_m))
path = trim(path) // "(" // trim(to_str(old_k)) // &
"," // trim(to_str(old_l)) // "," // &
trim(to_str(old_m)) // ")"
call find_offset(map, goal, next_univ, final, offset, path)
return
end if
end do
end do
end do
! ==================================================================
! HEXAGONAL LATTICES
type is (HexLattice)
! Write to the geometry stack
path = trim(path) // "->" // to_str(lat%id)
n_z = lat % n_axial
n_y = 2 * lat % n_rings - 1
n_x = 2 * lat % n_rings - 1
old_m = 1
old_l = 1
old_k = 1
! Loop over lattice coordinates
do m = 1, n_z
do l = 1, n_y
do k = 1, n_x
! This array position is never used
if (k + l < lat % n_rings + 1) then
cycle
! This array position is never used
else if (k + l > 3*lat % n_rings - 1) then
cycle
end if
if (final >= lat % offset(map, k, l, m) + offset) then
old_m = m
old_l = l
old_k = k
cycle
else
! Target is at this lattice position
lat_offset = lat % offset(map, old_k, old_l, old_m)
offset = offset + lat_offset
next_univ => universes(lat % universes(old_k, old_l, old_m))
path = trim(path) // "(" // &
trim(to_str(old_k - lat % n_rings)) // "," // &
trim(to_str(old_l - lat % n_rings)) // "," // &
trim(to_str(old_m)) // ")"
call find_offset(map, goal, next_univ, final, offset, path)
return
end if
end do
end do
end do
end select
end if
end if
end do
end subroutine find_offset
end module output

6
src/particle_header.F90
View file

@ -27,6 +27,9 @@ module particle_header
! Is this level rotated?
logical :: rotated = .false.
! Distributed Mapping Info
integer, pointer :: mapping(:) => null()
! Pointer to next (more local) set of coordinates
type(LocalCoord), pointer :: next => null()
end type LocalCoord
@ -102,6 +105,9 @@ contains
if (associated(coord % next)) call deallocate_coord(coord%next)
! deallocate original coordinate
if (associated(coord % mapping)) deallocate(coord % mapping)
! deallocate this coord
deallocate(coord)
end if

4
src/relaxng/tallies.rnc
View file

@ -23,9 +23,9 @@ element tallies {
attribute estimator { ( "analog" | "tracklength" ) })? &
element filter {
(element type { ( "cell" | "cellborn" | "material" | "universe" |
"surface" | "mesh" | "energy" | "energyout" ) } |
"surface" | "distribcell" | "mesh" | "energy" | "energyout" ) } |
attribute type { ( "cell" | "cellborn" | "material" | "universe" |
"surface" | "mesh" | "energy" | "energyout" ) }) &
"surface" | "distribcell" | "mesh" | "energy" | "energyout" ) }) &
(element bins { list { xsd:double+ } } |
attribute bins { list { xsd:double+ } })
}* &

2
src/relaxng/tallies.rng
View file

@ -147,6 +147,7 @@
<value>material</value>
<value>universe</value>
<value>surface</value>
<value>distribcell</value>
<value>mesh</value>
<value>energy</value>
<value>energyout</value>
@ -159,6 +160,7 @@
<value>material</value>
<value>universe</value>
<value>surface</value>
<value>distribcell</value>
<value>mesh</value>
<value>energy</value>
<value>energyout</value>

13
src/state_point.F90
View file

@ -249,6 +249,9 @@ contains
call sp % write_data(tally % filters(j) % type, "type", &
group="tallies/tally " // trim(to_str(tally % id)) // &
"/filter " // to_str(j))
call sp % write_data(tally % filters(j) % offset, "offset", &
group="tallies/tally " // trim(to_str(tally % id)) // &
"/filter " // to_str(j))
call sp % write_data(tally % filters(j) % n_bins, "n_bins", &
group="tallies/tally " // trim(to_str(tally % id)) // &
"/filter " // to_str(j))
@ -273,7 +276,11 @@ contains
! Set up nuclide bin array and then write
allocate(key_array(tally % n_nuclide_bins))
NUCLIDE_LOOP: do j = 1, tally % n_nuclide_bins
key_array(j) = tally % nuclide_bins(j)
if (tally % nuclide_bins(j) > 0) then
key_array(j) = nuclides(tally % nuclide_bins(j)) % zaid
else
key_array(j) = tally % nuclide_bins(j)
end if
end do NUCLIDE_LOOP
call sp % write_data(key_array, "nuclides", &
group="tallies/tally " // trim(to_str(tally % id)), &
@ -829,6 +836,9 @@ contains
call sp % read_data(tally % filters(j) % type, "type", &
group="tallies/tally " // trim(to_str(curr_key)) // &
"/filter " // to_str(j))
call sp % read_data(tally % filters(j) % offset, "offset", &
group="tallies/tally " // trim(to_str(curr_key)) // &
"/filter " // to_str(j))
call sp % read_data(tally % filters(j) % n_bins, "n_bins", &
group="tallies/tally " // trim(to_str(curr_key)) // &
"/filter " // to_str(j))
@ -863,6 +873,7 @@ contains
tally % nuclide_bins(j) = temp_array(j)
end if
end do NUCLIDE_LOOP
deallocate(temp_array)
! Write number of score bins, score bins, user score bins

21
src/tally.F90
View file

@ -3,6 +3,7 @@ module tally
use ace_header, only: Reaction
use constants
use error, only: fatal_error
use geometry_header
use global
use math, only: t_percentile, calc_pn, calc_rn
use mesh, only: get_mesh_bin, bin_to_mesh_indices, &
@ -10,7 +11,7 @@ module tally
mesh_intersects_2d, mesh_intersects_3d
use mesh_header, only: StructuredMesh
use output, only: header
use particle_header, only: LocalCoord, Particle
use particle_header, only: LocalCoord, Particle, deallocate_coord
use search, only: binary_search
use string, only: to_str
use tally_header, only: TallyResult, TallyMapItem, TallyMapElement
@ -1109,6 +1110,7 @@ contains
integer :: i ! loop index for filters
integer :: n ! number of bins for single filter
integer :: offset ! offset for distribcell
real(8) :: E ! particle energy
type(TallyObject), pointer :: t
type(StructuredMesh), pointer :: m
@ -1152,6 +1154,23 @@ contains
end do
nullify(coord)
case (FILTER_DISTRIBCELL)
! determine next distribcell bin
matching_bins(i) = NO_BIN_FOUND
coord => p % coord0
offset = 0
do while(associated(coord))
if (associated(coord % mapping)) then
offset = offset + coord % mapping(t % filters(i) % offset)
end if
if (t % filters(i) % int_bins(1) == coord % cell) then
matching_bins(i) = offset + 1
exit
end if
coord => coord % next
end do
nullify(coord)
case (FILTER_CELLBORN)
! determine next cellborn bin
matching_bins(i) = get_next_bin(FILTER_CELLBORN, &

1
src/tally_header.F90
View file

@ -54,6 +54,7 @@ module tally_header
type TallyFilter
integer :: type = NONE
integer :: n_bins = 0
integer :: offset = 0 ! Only used for distribcell filters
integer, allocatable :: int_bins(:)
real(8), allocatable :: real_bins(:) ! Only used for energy filters

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

@ -539,14 +539,14 @@ class CMFDFile(object):
self._write_matrices = write_matrices
def get_active_flush_subelement(self):
def create_active_flush_subelement(self):
if not self._active_flush is None:
element = ET.SubElement(self._cmfd_file, "active_flush")
element.text = '{0}'.format(str(self._active_flush))
def get_begin_subelement(self):
def create_begin_subelement(self):
if not self._begin is None:
element = ET.SubElement(self._cmfd_file, "begin")

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

@ -342,6 +342,10 @@ class StatePoint(object):
filter_type = self._get_int(
path='{0}{1}/type'.format(subbase, j))[0]
# Read the Filter offset
offset = self._get_int(
path='{0}{1}/offset'.format(subbase, j))[0]
n_bins = self._get_int(
path='{0}{1}/n_bins'.format(subbase, j))[0]
@ -366,6 +370,7 @@ class StatePoint(object):
# Create Filter object
filter = openmc.Filter(FILTER_TYPES[filter_type], bins)
filter.offset = offset
filter.stride = stride
filter.num_bins = n_bins
@ -684,6 +689,7 @@ class StatePoint(object):
nuclide_zaids = copy.deepcopy(tally._nuclides)
for nuclide_zaid in nuclide_zaids:
tally.remove_nuclide(nuclide_zaid)
if nuclide_zaid == -1:
tally.add_nuclide(openmc.Nuclide('total'))

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

@ -277,6 +277,12 @@ class Summary(object):
cell = openmc.Cell(cell_id=cell_id, name=name)
if fill_type == 'universe':
maps = self._f['geometry/cells'][key]['maps'][0]
if maps > 0:
offset = self._f['geometry/cells'][key]['offset'][...]
cell.set_offset(offset)
translated = self._f['geometry/cells'][key]['translated'][0]
if translated:
translation = \
@ -354,6 +360,11 @@ class Summary(object):
index = self._f['geometry/lattices'][key]['index'][0]
name = self._f['geometry/lattices'][key]['name'][0]
lattice_type = self._f['geometry/lattices'][key]['type'][...][0]
maps = self._f['geometry/lattices'][key]['maps'][0]
offset_size = self._f['geometry/lattices'][key]['offset_size'][0]
if offset_size > 0:
offsets = self._f['geometry/lattices'][key]['offsets'][...]
if lattice_type == 'rectangular':
dimension = self._f['geometry/lattices'][key]['dimension'][...]
@ -395,6 +406,11 @@ class Summary(object):
universes = universes[:,::-1,:]
lattice.universes = universes
if offset_size > 0:
offsets = np.swapaxes(offsets, 0, 1)
offsets = np.swapaxes(offsets, 1, 2)
lattice.offsets = offsets
# Add the Lattice to the global dictionary of all Lattices
self.lattices[index] = lattice
@ -430,6 +446,9 @@ class Summary(object):
universes[i,j,k] = \
self.get_universe_by_id(universe_ids[i,j,k])
if offset_size > 0:
lattice.offsets = offsets
# Add the Lattice to the global dictionary of all Lattices
self.lattices[index] = lattice
@ -514,15 +533,15 @@ class Summary(object):
subsubbase = '{0}/filter {1}'.format(subbase, j)
# Read filter type (e.g., "cell", "energy", etc.)
filter_type = self._f['{0}/type_name'.format(subsubbase)][0]
# Read filter type (e.g., "cell", "energy", etc.) integer code
filter_type = self._f['{0}/type'.format(subsubbase)][0]
# Read the filter bins
num_bins = self._f['{0}/n_bins'.format(subsubbase)][0]
bins = self._f['{0}/bins'.format(subsubbase)][...]
# Create Filter object
filter = openmc.Filter(filter_type, bins)
filter = openmc.Filter(openmc.FILTER_TYPES[filter_type], bins)
filter.num_bins = num_bins
# Add Filter to the Tally

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

@ -34,7 +34,7 @@ class Cell(object):
self._surfaces = {}
self._rotation = None
self._translation = None
self._offset = None
self._offsets = None
@property
@ -73,8 +73,8 @@ class Cell(object):
@property
def offset(self):
return self._offset
def offsets(self):
return self._offsets
@id.setter
@ -188,16 +188,16 @@ class Cell(object):
self._translation = translation
@offset.setter
def offset(self, offset):
@offsets.setter
def offsets(self, offsets):
if not isinstance(offset, (tuple, list, np.ndarray)):
msg = 'Unable to set offset {0} to Cell ID={1} since ' \
if not isinstance(offsets, (tuple, list, np.ndarray)):
msg = 'Unable to set offsets {0} to Cell ID={1} since ' \
'it is not a Python tuple/list or NumPy ' \
'array'.format(offset, self._id)
'array'.format(offsets, self._id)
raise ValueError(msg)
self._offset = offset
self._offsets = offsets
def add_surface(self, surface, halfspace):
@ -259,7 +259,7 @@ class Cell(object):
# If the Cell is filled by a Universe
elif self._type == 'fill':
offset = self._offset[filter_offset-1]
offset = self._offsets[filter_offset-1]
offset += self._fill.get_offset(path, filter_offset)
# If the Cell is filled by a Lattice
@ -268,12 +268,6 @@ class Cell(object):
return offset
# Make a recursive call to the Universe filling this Cell
offset = self._cells[cell_id].get_offset(path, filter_offset)
# Return the offset computed at all nested Universe levels
return offset
def get_all_nuclides(self):
@ -335,7 +329,7 @@ class Cell(object):
self._rotation)
string += '{0: <16}{1}{2}\n'.format('\tTranslation', '=\t',
self._translation)
string += '{0: <16}{1}{2}\n'.format('\tOffset', '=\t', self._offset)
string += '{0: <16}{1}{2}\n'.format('\tOffset', '=\t', self._offsets)
return string
@ -1170,7 +1164,7 @@ class HexLattice(Lattice):
raise ValueError(msg)
elif len(pitch) != 1 and len(pitch) != 2:
elif len(pitch) != 2 and len(pitch) != 3:
msg = 'Unable to set Lattice ID={0} pitch to {1} since it does ' \
'not contain 2 or 3 coordinates'.format(self._id, pitch)
raise ValueError(msg)

4
tests/run_tests.py
View file

@ -43,8 +43,8 @@ parser.add_option("-s", "--script", action="store_true", dest="script",
# Default compiler paths
FC='gfortran'
MPI_DIR='/opt/mpich/3.1.3-gnu'
HDF5_DIR='/opt/hdf5/1.8.14-gnu'
PHDF5_DIR='/opt/phdf5/1.8.14-gnu'
HDF5_DIR='/opt/hdf5/1.8.15-gnu'
PHDF5_DIR='/opt/phdf5/1.8.15-gnu'
# Script mode for extra capability
script_mode = False

24
tests/test_filter_distribcell/case-1/geometry.xml Normal file
View file

@ -0,0 +1,24 @@
<?xml version="1.0"?>
<geometry>
<cell id="1" fill="5" surfaces="1 -2 3 -4" />
<cell id="201" universe="2" material="1" surfaces="-5" />
<cell id="202" universe="2" material="2" surfaces="5" />
<lattice id="5">
<dimension>2 2</dimension>
<lower_left>-2.0 -2.0</lower_left>
<pitch>2.0 2.0</pitch>
<universes>
2 2
2 2
</universes>
</lattice>
<surface id="1" type="x-plane" coeffs="-2.0" boundary="vacuum" />
<surface id="2" type="x-plane" coeffs="2.0" boundary="vacuum" />
<surface id="3" type="y-plane" coeffs="-2.0" boundary="vacuum" />
<surface id="4" type="y-plane" coeffs="2.0" boundary="vacuum" />
<surface id="5" type="z-cylinder" coeffs="0.0 0.0 0.3" />
</geometry>

18
tests/test_filter_distribcell/case-1/materials.xml Normal file
View file

@ -0,0 +1,18 @@
<?xml version="1.0"?>
<materials>
<default_xs>71c</default_xs>
<!-- Definition of materials -->
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U-235" ao="1.0" />
</material>
<material id="2">
<density value="1.0" units="g/cc" />
<nuclide name="H-1" ao="2.0" />
<nuclide name="O-16" ao="1.0" />
</material>
</materials>

16
tests/test_filter_distribcell/case-1/settings.xml Normal file
View file

@ -0,0 +1,16 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>1</batches>
<inactive>0</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>-1 -1 -1 1 1 1</parameters>
</space>
</source>
</settings>

22
tests/test_filter_distribcell/case-1/tallies.xml Normal file
View file

@ -0,0 +1,22 @@
<?xml version="1.0"?>
<tallies>
<tally id="1">
<filter type="distribcell">
<bins>
201
</bins>
</filter>
<scores>total</scores>
</tally>
<tally id="2">
<filter type="cell">
<bins>
201
</bins>
</filter>
<scores>total</scores>
</tally>
</tallies>

30
tests/test_filter_distribcell/case-2/geometry.xml Normal file
View file

@ -0,0 +1,30 @@
<?xml version="1.0"?>
<geometry>
<cell id="1" fill="5" surfaces="1 -2 3 -4" />
<cell id="201" universe="21" material="1" surfaces="-5" />
<cell id="202" universe="21" material="2" surfaces="5" />
<cell id="203" universe="22" material="1" surfaces="-5" />
<cell id="204" universe="22" material="2" surfaces="5" />
<cell id="205" universe="23" material="1" surfaces="-5" />
<cell id="206" universe="23" material="2" surfaces="5" />
<cell id="207" universe="20" material="1" surfaces="-5" />
<cell id="208" universe="20" material="2" surfaces="5" />
<lattice id="5">
<dimension>2 2</dimension>
<lower_left>-2.0 -2.0</lower_left>
<pitch>2.0 2.0</pitch>
<universes>
20 21
22 23
</universes>
</lattice>
<surface id="1" type="x-plane" coeffs="-2.0" boundary="vacuum" />
<surface id="2" type="x-plane" coeffs="2.0" boundary="vacuum" />
<surface id="3" type="y-plane" coeffs="-2.0" boundary="vacuum" />
<surface id="4" type="y-plane" coeffs="2.0" boundary="vacuum" />
<surface id="5" type="z-cylinder" coeffs="0.0 0.0 0.3" />
</geometry>

18
tests/test_filter_distribcell/case-2/materials.xml Normal file
View file

@ -0,0 +1,18 @@
<?xml version="1.0"?>
<materials>
<default_xs>71c</default_xs>
<!-- Definition of materials -->
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U-235" ao="1.0" />
</material>
<material id="2">
<density value="1.0" units="g/cc" />
<nuclide name="H-1" ao="2.0" />
<nuclide name="O-16" ao="1.0" />
</material>
</materials>

16
tests/test_filter_distribcell/case-2/settings.xml Normal file
View file

@ -0,0 +1,16 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>1</batches>
<inactive>0</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>-1 -1 -1 1 1 1</parameters>
</space>
</source>
</settings>

13
tests/test_filter_distribcell/case-2/tallies.xml Normal file
View file

@ -0,0 +1,13 @@
<?xml version="1.0"?>
<tallies>
<tally id="1">
<filter type="cell">
<bins>
201 203 205 207
</bins>
</filter>
<scores>total</scores>
</tally>
</tallies>

181
tests/test_filter_distribcell/case-3/geometry.xml Normal file
View file

@ -0,0 +1,181 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="z-cylinder" coeffs="0. 0. 0.41" />
<surface id="2" type="z-cylinder" coeffs="0. 0. 0.475" />
<surface id="3" type="z-cylinder" coeffs="0. 0. 0.56" />
<surface id="4" type="z-cylinder" coeffs="0. 0. 0.62" />
<surface id="5" type="z-cylinder" coeffs="0. 0. 187.6" />
<surface id="6" type="z-cylinder" coeffs="0. 0. 209.0" />
<surface id="7" type="z-cylinder" coeffs="0. 0. 229.0" />
<surface id="8" type="z-cylinder" coeffs="0. 0. 249.0" boundary="vacuum" />
<surface id="31" type="z-plane" coeffs="-229.0" boundary="vacuum" />
<surface id="32" type="z-plane" coeffs="-199.0" />
<surface id="33" type="z-plane" coeffs="-193.0" />
<surface id="34" type="z-plane" coeffs="-183.0" />
<surface id="35" type="z-plane" coeffs="0.0" />
<surface id="36" type="z-plane" coeffs="183.0" />
<surface id="37" type="z-plane" coeffs="203.0" />
<surface id="38" type="z-plane" coeffs="215.0" />
<surface id="39" type="z-plane" coeffs="223.0" boundary="vacuum" />
<!-- All geometry on base universe -->
<cell id="1" fill="200" surfaces=" -6 34 -35" /> <!-- Lower core -->
<cell id="2" fill="201" surfaces=" -6 35 -36" /> <!-- Upper core -->
<cell id="3" material="8" surfaces=" -7 31 -32" /> <!-- Lower core plate region -->
<cell id="4" material="9" surfaces=" -5 32 -33" /> <!-- Bottom nozzle region -->
<cell id="5" material="12" surfaces=" -5 33 -34" /> <!-- Bottom FA region -->
<cell id="6" material="11" surfaces=" -5 36 -37" /> <!-- Top FA region -->
<cell id="7" material="10" surfaces=" -5 37 -38" /> <!-- Top nozzle region -->
<cell id="8" material="7" surfaces=" -7 38 -39" /> <!-- Upper plate region -->
<cell id="9" material="4" surfaces="6 -7 32 -38" /> <!-- Downcomer -->
<cell id="10" material="5" surfaces="7 -8 31 -39" /> <!-- RPV -->
<cell id="11" material="6" surfaces="5 -6 32 -34" /> <!-- Bottom of radial reflector -->
<cell id="12" material="7" surfaces="5 -6 36 -38" /> <!-- Top of radial reflector -->
<!-- Fuel pin, cladding, cold water -->
<cell id="21" universe="1" material="1" surfaces="-1" />
<cell id="22" universe="1" material="2" surfaces="1 -2" />
<cell id="23" universe="1" material="3" surfaces="2" />
<!-- Instrumentation guide tube -->
<cell id="24" universe="2" material="3" surfaces="-3" />
<cell id="25" universe="2" material="2" surfaces="3 -4" />
<cell id="26" universe="2" material="3" surfaces="4" />
<!-- Fuel pin, cladding, hot water -->
<cell id="27" universe="3" material="1" surfaces="-1" />
<cell id="28" universe="3" material="2" surfaces="1 -2" />
<cell id="29" universe="3" material="4" surfaces="2" />
<!-- Instrumentation guide tube -->
<cell id="30" universe="4" material="4" surfaces="-3" />
<cell id="31" universe="4" material="2" surfaces="3 -4" />
<cell id="32" universe="4" material="4" surfaces="4" />
<!-- cell for water assembly (cold) -->
<cell id="50" universe="5" material="4" surfaces="34 -35" />
<!-- containing cell for fuel assembly -->
<cell id="60" universe="6" fill="100" surfaces="34 -35" />
<!-- cell for water assembly (hot) -->
<cell id="70" universe="7" material="3" surfaces="35 -36" />
<!-- containing cell for fuel assembly -->
<cell id="80" universe="8" fill="101" surfaces="35 -36" />
<!-- Fuel Assembly (Lower Half) -->
<lattice id="100">
<dimension>17 17</dimension>
<lower_left>-10.71 -10.71</lower_left>
<pitch>1.26 1.26</pitch>
<universes>
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1
1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1
1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
</universes>
</lattice>
<!-- Fuel Assembly (Upper Half) -->
<lattice id="101">
<dimension>17 17</dimension>
<lower_left>-10.71 -10.71</lower_left>
<pitch>1.26 1.26</pitch>
<universes>
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3
3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3
3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
</universes>
</lattice>
<!-- Core Lattice (Lower Half) -->
<lattice id="200">
<dimension>21 21</dimension>
<lower_left>-224.91 -224.91</lower_left>
<pitch>21.42 21.42</pitch>
<universes>
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5
5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5
5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5
5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5
5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5
5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
</universes>
</lattice>
<!-- Core Lattice (Upper Half) -->
<lattice id="201">
<dimension>21 21</dimension>
<lower_left>-224.91 -224.91</lower_left>
<pitch>21.42 21.42</pitch>
<universes>
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7
7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7
7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7
7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7
7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7
7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
</universes>
</lattice>
</geometry>

125
tests/test_filter_distribcell/case-3/materials.xml Normal file
View file

@ -0,0 +1,125 @@
<?xml version="1.0"?>
<materials>
<default_xs>71c</default_xs>
<!-- Fuel composition -->
<material id="1">
<density value="10.062" units="g/cm3" />
<nuclide name="U-235" ao="4.8218e-4" />
<nuclide name="U-238" ao="2.1504e-2" />
<nuclide name="O-16" ao="4.5737e-2" />
</material>
<!-- Cladding composition -->
<material id="2">
<density value="5.77" units="g/cm3" />
<nuclide name="Zr-90" ao="0.5145" />
<nuclide name="Zr-91" ao="0.1122" />
<nuclide name="Zr-92" ao="0.1715" />
<nuclide name="Zr-94" ao="0.1738" />
<nuclide name="Zr-96" ao="0.0280" />
</material>
<!-- Cold borated water -->
<material id="3">
<density value="0.07416" units="atom/b-cm" />
<nuclide name="H-1" ao="2.0" />
<nuclide name="O-16" ao="1.0" />
<nuclide name="B-10" ao="6.490e-4" />
<nuclide name="B-11" ao="2.689e-3" />
</material>
<!-- Hot borated water -->
<material id="4">
<density value="0.06614" units="atom/b-cm" />
<nuclide name="H-1" ao="2.0" />
<nuclide name="O-16" ao="1.0" />
<nuclide name="B-10" ao="6.490e-4" />
<nuclide name="B-11" ao="2.689e-3" />
</material>
<!-- RPV Composition -->
<material id="5">
<density value="7.9" units="g/cm3" />
<nuclide name="Zr-90" ao="0.5145" />
<nuclide name="Zr-91" ao="0.1122" />
<nuclide name="Zr-92" ao="0.1715" />
<nuclide name="Zr-94" ao="0.1738" />
<nuclide name="Zr-96" ao="0.0280" />
</material>
<!-- Lower radial reflector -->
<material id="6">
<density value="4.32" units="g/cm3" />
<nuclide name="H-1" wo="0.0095661" />
<nuclide name="O-16" wo="0.0759107" />
<nuclide name="B-10" wo="3.08409e-5" />
<nuclide name="B-11" wo="1.40499e-4" />
</material>
<!-- Upper radial reflector / Top plate region -->
<material id="7">
<density value="4.28" units="g/cm3" />
<nuclide name="H-1" wo="0.0086117" />
<nuclide name="O-16" wo="0.0683369" />
<nuclide name="B-10" wo="2.77638e-5" />
<nuclide name="B-11" wo="1.26481e-4" />
</material>
<!-- Bottom plate region -->
<material id="8">
<density value="7.184" units="g/cm3" />
<nuclide name="H-1" wo="0.0011505" />
<nuclide name="O-16" wo="0.0091296" />
<nuclide name="B-10" wo="3.70915e-6" />
<nuclide name="B-11" wo="1.68974e-5" />
</material>
<!-- Bottom nozzle region -->
<material id="9">
<density value="2.53" units="g/cm3" />
<nuclide name="H-1" wo="0.0245014" />
<nuclide name="O-16" wo="0.1944274" />
<nuclide name="B-10" wo="7.89917e-5" />
<nuclide name="B-11" wo="3.59854e-4" />
</material>
<!-- Top nozzle region -->
<material id="10">
<density value="1.746" units="g/cm3" />
<nuclide name="H-1" wo="0.0358870" />
<nuclide name="O-16" wo="0.2847761" />
<nuclide name="B-10" wo="1.15699e-4" />
<nuclide name="B-11" wo="5.27075e-4" />
</material>
<!-- Top of Fuel Assemblies -->
<material id="11">
<density value="3.044" units="g/cm3" />
<nuclide name="H-1" wo="0.0162913" />
<nuclide name="O-16" wo="0.1292776" />
<nuclide name="B-10" wo="5.25228e-5" />
<nuclide name="B-11" wo="2.39272e-4" />
<nuclide name="Zr-90" wo="0.43313403903" />
<nuclide name="Zr-91" wo="0.09549277374" />
<nuclide name="Zr-92" wo="0.14759527104" />
<nuclide name="Zr-94" wo="0.15280552077" />
<nuclide name="Zr-96" wo="0.02511169542" />
</material>
<!-- Bottom of Fuel Assemblies -->
<material id="12">
<density value="1.762" units="g/cm3" />
<nuclide name="H-1" wo="0.0292856" />
<nuclide name="O-16" wo="0.2323919" />
<nuclide name="B-10" wo="9.44159e-5" />
<nuclide name="B-11" wo="4.30120e-4" />
<nuclide name="Zr-90" wo="0.3741373658" />
<nuclide name="Zr-91" wo="0.0824858164" />
<nuclide name="Zr-92" wo="0.1274914944" />
<nuclide name="Zr-94" wo="0.1319920622" />
<nuclide name="Zr-96" wo="0.0216912612" />
</material>
</materials>

21
tests/test_filter_distribcell/case-3/settings.xml Normal file
View file

@ -0,0 +1,21 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>3</batches>
<inactive>0</inactive>
<particles>100</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>
-160 -160 -183
160 160 183
</parameters>
</space>
</source>
<output tallies="false" />
</settings>

58
tests/test_filter_distribcell/case-3/tallies.xml Normal file
View file

@ -0,0 +1,58 @@
<?xml version="1.0"?>
<tallies>
<tally id="1">
<filter type="cell">
<bins>
1
</bins>
</filter>
<scores>total</scores>
</tally>
<tally id="2">
<filter type="distribcell">
<bins>
1
</bins>
</filter>
<scores>total</scores>
</tally>
<tally id="3">
<filter type="cell">
<bins>
60
</bins>
</filter>
<scores>total</scores>
</tally>
<tally id="4">
<filter type="distribcell">
<bins>
60
</bins>
</filter>
<scores>total</scores>
</tally>
<tally id="5">
<filter type="cell">
<bins>
27
</bins>
</filter>
<scores>total</scores>
</tally>
<tally id="6">
<filter type="distribcell">
<bins>
27
</bins>
</filter>
<scores>total</scores>
</tally>
</tallies>

111
tests/test_filter_distribcell/results.py Normal file
View file

@ -0,0 +1,111 @@
#!/usr/bin/env python
import sys
import numpy as np
# import statepoint
sys.path.insert(0, '../../src/utils')
from openmc.statepoint import StatePoint
from openmc import Filter
# read in statepoint file
sp3 = StatePoint(sys.argv[1])
sp2 = StatePoint(sys.argv[2])
sp1 = StatePoint(sys.argv[3])
sp1.read_results()
sp2.read_results()
sp3.read_results()
sp1.compute_stdev()
sp2.compute_stdev()
sp3.compute_stdev()
# analyze sp1
# compare distrib sum to cell
sp1_t1 = sp1._tallies[1]
sp1_t2 = sp1._tallies[2]
distribcell_filter = sp1_t1.find_filter(filter_type='distribcell', bins=[2])
sp1_t1_d1 = sp1_t1.get_value(score='total', filters=[distribcell_filter],
filter_bins=[0], value='mean')
sp1_t1_d2 = sp1_t1.get_value(score='total', filters=[distribcell_filter],
filter_bins=[1], value='mean')
sp1_t1_d3 = sp1_t1.get_value(score='total', filters=[distribcell_filter],
filter_bins=[2], value='mean')
sp1_t1_d4 = sp1_t1.get_value(score='total', filters=[distribcell_filter],
filter_bins=[3], value='mean')
sp1_t1_sum = sp1_t1_d1 + sp1_t1_d2 + sp1_t1_d3 + sp1_t1_d4
cell_filter = sp1_t2.find_filter(filter_type='cell', bins=[2])
sp1_t2_c201 = sp1_t2.get_value(score='total', filters=[cell_filter],
filter_bins=[2], value='mean')
# analyze sp2
sp2_t1 = sp2._tallies[1]
cell_filter = sp2_t1.find_filter(filter_type='cell', bins=[2,4,6,8])
sp2_t1_c201 = sp2_t1.get_value(score='total', filters=[cell_filter],
filter_bins=[2], value='mean')
sp2_t1_c203 = sp2_t1.get_value(score='total', filters=[cell_filter],
filter_bins=[4], value='mean')
sp2_t1_c205 = sp2_t1.get_value(score='total', filters=[cell_filter],
filter_bins=[6], value='mean')
sp2_t1_c207 = sp2_t1.get_value(score='total', filters=[cell_filter],
filter_bins=[8], value='mean')
# analyze sp3
sp3_t1 = sp3._tallies[1]
sp3_t2 = sp3._tallies[2]
sp3_t3 = sp3._tallies[3]
sp3_t4 = sp3._tallies[4]
sp3_t5 = sp3._tallies[5]
sp3_t6 = sp3._tallies[6]
cell_filter = sp3_t1.find_filter(filter_type='cell', bins=[1])
distribcell_filter = sp3_t2.find_filter(filter_type='distribcell', bins=[1])
sp3_t1_c1 = sp3_t1.get_value(score='total', filters=[cell_filter],
filter_bins=[1], value='mean')
sp3_t2_d1 = sp3_t2.get_value(score='total', filters=[distribcell_filter],
filter_bins=[0], value='mean')
cell_filter = sp3_t3.find_filter(filter_type='cell', bins=[26])
distribcell_filter = sp3_t4.find_filter(filter_type='distribcell', bins=[26])
sp3_t3_c60 = sp3_t3.get_value(score='total', filters=[cell_filter],
filter_bins=[26], value='mean')
sp3_t4_d = 0
for i in range(241):
sp3_t4_d += sp3_t4.get_value(score='total', filters=[distribcell_filter],
filter_bins=[i], value='mean')
cell_filter = sp3_t5.find_filter(filter_type='cell', bins=[19])
distribcell_filter = sp3_t6.find_filter(filter_type='distribcell', bins=[19])
sp3_t5_c27 = sp3_t5.get_value(score='total', filters=[cell_filter],
filter_bins=[19], value='mean')
sp3_t6_d = 0
for i in range(63624):
sp3_t6_d += sp3_t6.get_value(score='total', filters=[distribcell_filter],
filter_bins=[i], value='mean')
# set up output string
outstr = ''
outstr += "{0:12.6E}\n".format(sp1_t1_sum)
outstr += "{0:12.6E}\n".format(sp1_t2_c201)
outstr += "{0:12.6E}\n".format(sp2_t1_c201)
outstr += "{0:12.6E}\n".format(sp1_t1_d4)
outstr += "{0:12.6E}\n".format(sp2_t1_c203)
outstr += "{0:12.6E}\n".format(sp1_t1_d1)
outstr += "{0:12.6E}\n".format(sp2_t1_c205)
outstr += "{0:12.6E}\n".format(sp1_t1_d3)
outstr += "{0:12.6E}\n".format(sp2_t1_c207)
outstr += "{0:12.6E}\n".format(sp1_t1_d2)
outstr += "{0:12.6E}\n".format(sp3_t1_c1)
outstr += "{0:12.6E}\n".format(sp3_t2_d1)
outstr += "{0:12.6E}\n".format(sp3_t3_c60)
outstr += "{0:12.6E}\n".format(sp3_t4_d)
outstr += "{0:12.6E}\n".format(sp3_t5_c27)
outstr += "{0:12.6E}\n".format(sp3_t6_d)
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

16
tests/test_filter_distribcell/results_true.dat Normal file
View file

@ -0,0 +1,16 @@
3.296526E-02
3.296526E-02
7.326285E-03
7.326285E-03
8.565980E-03
8.565980E-03
9.027116E-03
9.027116E-03
8.045879E-03
8.045879E-03
1.820738E+01
1.820738E+01
1.770306E+01
1.770306E+01
2.109466E+00
2.109466E+00

108
tests/test_filter_distribcell/test_filter_distribcell.py Normal file
View file

@ -0,0 +1,108 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
import glob
from optparse import OptionParser
from shutil import copyfile
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()
def test_run():
if opts.mpi_exec != '':
opts.mpi_exe = os.path.abspath(opts.mpi_exec)
opts.exe = os.path.abspath(opts.exe)
# Call 1
os.chdir('case-1')
cwd = os.getcwd()
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
# Call 2
os.chdir('../case-2')
cwd = os.getcwd()
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
# Call 3
os.chdir('../case-3')
cwd = os.getcwd()
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
os.chdir('..')
def test_created_statepoint():
cwd = os.getcwd()
statepoint1 = glob.glob(cwd + '/case-1/statepoint.1.*')
statepoint2 = glob.glob(cwd + '/case-2/statepoint.1.*')
statepoint3 = glob.glob(cwd + '/case-3/statepoint.3.*')
assert len(statepoint1) == 1
assert len(statepoint2) == 1
assert len(statepoint3) == 1
string1 = statepoint1.pop()
string2 = statepoint2.pop()
string3 = statepoint3.pop()
assert string1.endswith('binary') or string1.endswith('h5')
assert string2.endswith('binary') or string2.endswith('h5')
assert string3.endswith('binary') or string3.endswith('h5')
def test_results():
cwd = os.getcwd()
statepoint = list()
statepoint.append(glob.glob(cwd + '/case-1/statepoint.1.*'))
statepoint.append(glob.glob(cwd + '/case-2/statepoint.1.*'))
statepoint.append(glob.glob(cwd + '/case-3/statepoint.3.*'))
call(['python', 'results.py', statepoint.pop()[0], statepoint.pop()[0], statepoint.pop()[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare
def teardown():
cwd = os.getcwd()
output = glob.glob(cwd + '/statepoint.*')
output.append(cwd + '/results_test.dat')
for f in output:
if os.path.exists(str(f)):
os.remove(str(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()

10
tests/travis_install.sh
View file

@ -15,9 +15,9 @@ if [ "$TRAVIS_PULL_REQUEST" != "false" ]; then
cd ..
# Build PHDF5
wget -q http://www.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8.14/src/hdf5-1.8.14.tar.gz
tar -xzvf hdf5-1.8.14.tar.gz >/dev/null 2>&1
mv hdf5-1.8.14 phdf5-1.8.14; cd phdf5-1.8.14
wget -q http://www.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8.15/src/hdf5-1.8.15.tar.gz
tar -xzvf hdf5-1.8.15.tar.gz >/dev/null 2>&1
mv hdf5-1.8.15 phdf5-1.8.15; cd phdf5-1.8.15
CC=$PWD/../mpich_install/bin/mpicc FC=$PWD/../mpich_install/bin/mpif90 \
./configure \
--prefix=$PWD/../phdf5_install -q --enable-fortran \
@ -27,8 +27,8 @@ if [ "$TRAVIS_PULL_REQUEST" != "false" ]; then
cd ..
# Build HDF5
tar -xzvf hdf5-1.8.14.tar.gz >/dev/null 2>&1
cd hdf5-1.8.14
tar -xzvf hdf5-1.8.15.tar.gz >/dev/null 2>&1
cd hdf5-1.8.15
CC=gcc FC=gfortran ./configure --prefix=$PWD/../hdf5_install -q \
--enable-fortran --enable-fortran2003
make -j >/dev/null 2>&1

17
tests/update_results.py
View file

@ -56,9 +56,8 @@ for adir in sorted(folders):
for i in range(35 - sz):
print('.', end="")
# Handle source file test separately since it requires running OpenMC twice
if adir == 'test_source_file':
# Handle source file test separately since it requires running OpenMC
# twice
if os.path.exists('results_error.dat'):
os.remove('results_error.dat')
proc = Popen(['python', 'test_source_file.py', '--exe', openmc_exe],
@ -70,6 +69,20 @@ for adir in sorted(folders):
os.chdir('..')
continue
# Handle distribcell test separately since it requires running OpenMC 3x
elif adir == 'test_filter_distribcell':
if os.path.exists('results_error.dat'):
os.remove('results_error.dat')
proc = Popen(['python', 'test_filter_distribcell.py',
'--exe', openmc_exe], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()
if os.path.exists('results_error.dat'):
print('renamed results_error!!!')
os.rename('results_error.dat', 'results_true.dat')
print(BOLD + OKGREEN + "[OK]" + ENDC)
os.chdir('..')
continue
# Run openmc
proc = Popen([openmc_exe], stderr=STDOUT, stdout=PIPE)
returncode = proc.wait()