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

Merge pull request #678 from smharper/oo_filters

Browse source 
Object-oriented tally filters
This commit is contained in:
Paul Romano 2016-08-02 10:16:10 -05:00 • committed by GitHub
commit ab432ac151
24 changed files with 2816 additions and 2342 deletions
Download patch file Download diff file Expand all files Collapse all files

5
openmc/filter.py
View file

@ -162,6 +162,11 @@ class Filter(object):
if not isinstance(bins, Iterable):
bins = [bins]
# If the bin is 0D numpy array, promote to 1D
elif isinstance(bins, np.ndarray):
if bins.shape == ():
bins.shape = (1,)
# If the bins are in a collection, convert it to a list
else:
bins = list(bins)

59
src/cmfd_data.F90
View file

@ -6,6 +6,7 @@ module cmfd_data
!==============================================================================
use constants
use tally_filter, only: MeshFilter
implicit none
private
@ -94,7 +95,10 @@ contains
! Associate tallies and mesh
t => cmfd_tallies(1)
i_mesh = t % filters(t % find_filter(FILTER_MESH)) % int_bins(1)
select type(filt => t % filters(t % find_filter(FILTER_MESH)) % obj)
type is (MeshFilter)
i_mesh = filt % mesh
end select
m => meshes(i_mesh)
! Set mesh widths
@ -109,7 +113,10 @@ contains
! Associate tallies and mesh
t => cmfd_tallies(ital)
i_mesh = t % filters(t % find_filter(FILTER_MESH)) % int_bins(1)
select type(filt => t % filters(t % find_filter(FILTER_MESH)) % obj)
type is (MeshFilter)
i_mesh = filt % mesh
end select
m => meshes(i_mesh)
i_filter_mesh = t % find_filter(FILTER_MESH)
@ -138,7 +145,7 @@ contains
TALLY: if (ital == 1) then
! Reset all bins to 1
matching_bins(1:t%n_filters) = 1
matching_bins(1:size(t % filters)) = 1
! Set ijk as mesh indices
ijk = (/ i, j, k /)
@ -152,7 +159,8 @@ contains
end if
! Calculate score index from bins
score_index = sum((matching_bins(1:t%n_filters) - 1) * t%stride) + 1
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t%stride) + 1
! Get flux
flux = t % results(1,score_index) % sum
@ -181,7 +189,7 @@ contains
INGROUP: do g = 1, ng
! Reset all bins to 1
matching_bins(1:t%n_filters) = 1
matching_bins(1:size(t % filters)) = 1
! Set ijk as mesh indices
ijk = (/ i, j, k /)
@ -198,7 +206,8 @@ contains
end if
! Calculate score index from bins
score_index = sum((matching_bins(1:t%n_filters) - 1) * t%stride) + 1
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t%stride) + 1
! Get scattering
cmfd % scattxs(h,g,i,j,k) = t % results(1,score_index) % sum /&
@ -220,7 +229,7 @@ contains
else if (ital == 3) then
! Initialize and filter for energy
matching_bins(1:t%n_filters) = 1
matching_bins(1:size(t % filters)) = 1
if (i_filter_ein > 0) then
matching_bins(i_filter_ein) = ng - h + 1
end if
@ -229,60 +238,72 @@ contains
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i-1, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_RIGHT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t%stride) + 1 ! outgoing
cmfd % current(1,h,i,j,k) = t % results(1,score_index) % sum
matching_bins(i_filter_surf) = OUT_RIGHT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! incoming
cmfd % current(2,h,i,j,k) = t % results(1,score_index) % sum
! Right surface
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_RIGHT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! incoming
cmfd % current(3,h,i,j,k) = t % results(1,score_index) % sum
matching_bins(i_filter_surf) = OUT_RIGHT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! outgoing
cmfd % current(4,h,i,j,k) = t % results(1,score_index) % sum
! Back surface
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j-1, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_FRONT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! outgoing
cmfd % current(5,h,i,j,k) = t % results(1,score_index) % sum
matching_bins(i_filter_surf) = OUT_FRONT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! incoming
cmfd % current(6,h,i,j,k) = t % results(1,score_index) % sum
! Front surface
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_FRONT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! incoming
cmfd % current(7,h,i,j,k) = t % results(1,score_index) % sum
matching_bins(i_filter_surf) = OUT_FRONT
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! outgoing
cmfd % current(8,h,i,j,k) = t % results(1,score_index) % sum
! Bottom surface
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k-1 /) + 1, .true.)
matching_bins(i_filter_surf) = IN_TOP
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! outgoing
cmfd % current(9,h,i,j,k) = t % results(1,score_index) % sum
matching_bins(i_filter_surf) = OUT_TOP
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! incoming
cmfd % current(10,h,i,j,k) = t % results(1,score_index) % sum
! Top surface
matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, &
(/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_TOP
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! incoming
cmfd % current(11,h,i,j,k) = t % results(1,score_index) % sum
matching_bins(i_filter_surf) = OUT_TOP
score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing
score_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1 ! outgoing
cmfd % current(12,h,i,j,k) = t % results(1,score_index) % sum
end if TALLY

97
src/cmfd_input.F90
View file

@ -267,14 +267,16 @@ contains
use mesh_header, only: RegularMesh
use string
use tally, only: setup_active_cmfdtallies
use tally_header, only: TallyObject, TallyFilter
use tally_header, only: TallyObject
use tally_filter_header
use tally_filter
use tally_initialize, only: add_tallies
use xml_interface
type(Node), pointer :: doc ! pointer to XML doc info
character(MAX_LINE_LEN) :: temp_str ! temp string
integer :: i ! loop counter
integer :: i, j ! loop counter
integer :: n ! size of arrays in mesh specification
integer :: ng ! number of energy groups (default 1)
integer :: n_filters ! number of filters
@ -283,7 +285,7 @@ contains
real(8) :: rarray3(3) ! temp double array
type(TallyObject), pointer :: t
type(RegularMesh), pointer :: m
type(TallyFilter) :: filters(N_FILTER_TYPES) ! temporary filters
type(TallyFilterContainer) :: filters(N_FILTER_TYPES) ! temporary filters
type(Node), pointer :: node_mesh
! Set global variables if they are 0 (this can happen if there is no tally
@ -410,21 +412,25 @@ contains
! Set up mesh filter
n_filters = 1
filters(n_filters) % type = FILTER_MESH
filters(n_filters) % n_bins = product(m % dimension)
allocate(filters(n_filters) % int_bins(1))
filters(n_filters) % int_bins(1) = n_user_meshes + 1
allocate(MeshFilter :: filters(n_filters) % obj)
select type (filt => filters(n_filters) % obj)
type is (MeshFilter)
filt % n_bins = product(m % dimension)
filt % mesh = n_user_meshes + 1
end select
t % find_filter(FILTER_MESH) = n_filters
! Read and set incoming energy mesh filter
if (check_for_node(node_mesh, "energy")) then
n_filters = n_filters + 1
filters(n_filters) % type = FILTER_ENERGYIN
ng = get_arraysize_double(node_mesh, "energy")
filters(n_filters) % n_bins = ng - 1
allocate(filters(n_filters) % real_bins(ng))
call get_node_array(node_mesh, "energy", &
filters(n_filters) % real_bins)
allocate(EnergyFilter :: filters(n_filters) % obj)
select type (filt => filters(n_filters) % obj)
type is (EnergyFilter)
ng = get_arraysize_double(node_mesh, "energy")
filt % n_bins = ng - 1
allocate(filt % bins(ng))
call get_node_array(node_mesh, "energy", filt % bins)
end select
t % find_filter(FILTER_ENERGYIN) = n_filters
end if
@ -448,9 +454,10 @@ contains
t % type = TALLY_VOLUME
! Allocate and set filters
t % n_filters = n_filters
allocate(t % filters(n_filters))
t % filters = filters(1:n_filters)
do j = 1, n_filters
call move_alloc(filters(j) % obj, t % filters(j) % obj)
end do
! Allocate scoring bins
allocate(t % score_bins(3))
@ -481,23 +488,22 @@ contains
! read and set outgoing energy mesh filter
if (check_for_node(node_mesh, "energy")) then
n_filters = n_filters + 1
filters(n_filters) % type = FILTER_ENERGYOUT
ng = get_arraysize_double(node_mesh, "energy")
filters(n_filters) % n_bins = ng - 1
allocate(filters(n_filters) % real_bins(ng))
call get_node_array(node_mesh, "energy", &
filters(n_filters) % real_bins)
allocate(EnergyoutFilter :: filters(n_filters) % obj)
select type (filt => filters(n_filters) % obj)
type is (EnergyoutFilter)
ng = get_arraysize_double(node_mesh, "energy")
filt % n_bins = ng - 1
allocate(filt % bins(ng))
call get_node_array(node_mesh, "energy", filt % bins)
end select
t % find_filter(FILTER_ENERGYOUT) = n_filters
end if
! Allocate and set filters
t % n_filters = n_filters
allocate(t % filters(n_filters))
t % filters = filters(1:n_filters)
! deallocate filters bins array
if (check_for_node(node_mesh, "energy")) &
deallocate(filters(n_filters) % real_bins)
do j = 1, n_filters
call move_alloc(filters(j) % obj, t % filters(j) % obj)
end do
! Allocate macro reactions
allocate(t % score_bins(2))
@ -522,25 +528,25 @@ contains
! Add extra filter for surface
n_filters = n_filters + 1
filters(n_filters) % type = FILTER_SURFACE
filters(n_filters) % n_bins = 2 * m % n_dimension
allocate(filters(n_filters) % int_bins(2 * m % n_dimension))
if (m % n_dimension == 2) then
filters(n_filters) % int_bins = (/ IN_RIGHT, OUT_RIGHT, IN_FRONT, &
OUT_FRONT /)
elseif (m % n_dimension == 3) then
filters(n_filters) % int_bins = (/ IN_RIGHT, OUT_RIGHT, IN_FRONT, &
OUT_FRONT, IN_TOP, OUT_TOP /)
end if
allocate(SurfaceFilter :: filters(n_filters) % obj)
select type(filt => filters(n_filters) % obj)
type is(SurfaceFilter)
filt % n_bins = 2 * m % n_dimension
allocate(filt % surfaces(2 * m % n_dimension))
if (m % n_dimension == 2) then
filt % surfaces = (/ IN_RIGHT, OUT_RIGHT, IN_FRONT, OUT_FRONT /)
elseif (m % n_dimension == 3) then
filt % surfaces = (/ IN_RIGHT, OUT_RIGHT, IN_FRONT, OUT_FRONT, &
IN_TOP, OUT_TOP /)
end if
end select
t % find_filter(FILTER_SURFACE) = n_filters
! Allocate and set filters
t % n_filters = n_filters
allocate(t % filters(n_filters))
t % filters = filters(1:n_filters)
! Deallocate filters bins array
deallocate(filters(n_filters) % int_bins)
do j = 1, n_filters
call move_alloc(filters(j) % obj, t % filters(j) % obj)
end do
! Allocate macro reactions
allocate(t % score_bins(1))
@ -558,15 +564,10 @@ contains
! We need to increase the dimension by one since we also need
! currents coming into and out of the boundary mesh cells.
i_filter_mesh = t % find_filter(FILTER_MESH)
t % filters(i_filter_mesh) % n_bins = product(m % dimension + 1)
t % filters(i_filter_mesh) % obj % n_bins = product(m % dimension + 1)
end if
! Deallocate filter bins
deallocate(filters(1) % int_bins)
if (check_for_node(node_mesh, "energy")) &
deallocate(filters(2) % real_bins)
end do
! Put cmfd tallies into active tally array and turn tallies on

4
src/geometry.F90
View file

@ -499,7 +499,7 @@ contains
end if
! Set previous coordinate going slightly past surface crossing
p % last_xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
! Diagnostic message
if (verbosity >= 10 .or. trace) then
@ -563,7 +563,7 @@ contains
end if
! Set previous coordinate going slightly past surface crossing
p % last_xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
! Diagnostic message
if (verbosity >= 10 .or. trace) then

11
src/global.F90
View file

@ -14,7 +14,7 @@ module global
use set_header, only: SetInt
use surface_header, only: SurfaceContainer
use source_header, only: SourceDistribution
use tally_header, only: TallyObject, TallyMap, TallyResult
use tally_header, only: TallyObject, TallyResult
use trigger_header, only: KTrigger
use timer_header, only: Timer
@ -137,6 +137,7 @@ module global
type(RegularMesh), allocatable, target :: meshes(:)
type(TallyObject), allocatable, target :: tallies(:)
integer, allocatable :: matching_bins(:)
real(8), allocatable :: filter_weights(:)
! Pointers for different tallies
type(TallyObject), pointer :: user_tallies(:) => null()
@ -176,9 +177,6 @@ module global
!$omp threadprivate(global_tally_collision, global_tally_absorption, &
!$omp& global_tally_tracklength, global_tally_leakage)
! Tally map structure
type(TallyMap), allocatable :: tally_maps(:)
integer :: n_meshes = 0 ! # of structured meshes
integer :: n_user_meshes = 0 ! # of structured user meshes
integer :: n_tallies = 0 ! # of tallies
@ -443,7 +441,8 @@ module global
type(Nuclide0K), allocatable, target :: nuclides_0K(:) ! 0K nuclides info
!$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, &
!$omp& trace, thread_id, current_work, matching_bins)
!$omp& trace, thread_id, current_work, matching_bins, &
!$omp& filter_weights)
contains
@ -509,7 +508,7 @@ contains
if (allocated(meshes)) deallocate(meshes)
if (allocated(tallies)) deallocate(tallies)
if (allocated(matching_bins)) deallocate(matching_bins)
if (allocated(tally_maps)) deallocate(tally_maps)
if (allocated(filter_weights)) deallocate(filter_weights)
! Deallocate fission and source bank and entropy
!$omp parallel

111
src/initialize.F90
View file

@ -22,8 +22,9 @@ module initialize
use state_point, only: load_state_point
use string, only: to_str, starts_with, ends_with, str_to_int
use summary, only: write_summary
use tally_header, only: TallyObject, TallyResult, TallyFilter
use tally_header, only: TallyObject, TallyResult
use tally_initialize,only: configure_tallies
use tally_filter
use tally, only: init_tally_routines
#ifdef MPI
@ -711,76 +712,15 @@ contains
! =======================================================================
! ADJUST INDICES FOR EACH TALLY FILTER
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
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_SURFACE)
FILTER_LOOP: do j = 1, size(t % filters)
select type(filt => t % filters(j) % obj)
type is (SurfaceFilter)
! Check if this is a surface filter only for surface currents
if (any(t%score_bins == SCORE_CURRENT)) cycle FILTER_LOOP
do k = 1, t%filters(j)%n_bins
id = t%filters(j)%int_bins(k)
if (surface_dict%has_key(id)) then
t%filters(j)%int_bins(k) = surface_dict%get_key(id)
else
call fatal_error("Could not find surface " // trim(to_str(id)) &
&// " specified on tally " // trim(to_str(t%id)))
end if
end do
case (FILTER_UNIVERSE)
do k = 1, t%filters(j)%n_bins
id = t%filters(j)%int_bins(k)
if (universe_dict%has_key(id)) then
t%filters(j)%int_bins(k) = universe_dict%get_key(id)
else
call fatal_error("Could not find universe " // trim(to_str(id)) &
&// " specified on tally " // trim(to_str(t%id)))
end if
end do
case (FILTER_MATERIAL)
do k = 1, t%filters(j)%n_bins
id = t%filters(j)%int_bins(k)
if (material_dict%has_key(id)) then
t%filters(j)%int_bins(k) = material_dict%get_key(id)
else
call fatal_error("Could not find material " // trim(to_str(id)) &
&// " specified on tally " // trim(to_str(t%id)))
end if
end do
case (FILTER_MESH)
! The mesh filter already has been set to the index in meshes rather
! than the user-specified id, so it doesn't need to be changed.
if (.not. any(t % score_bins == SCORE_CURRENT)) &
call filt % initialize()
class default
call filt % initialize()
end select
end do FILTER_LOOP
@ -894,14 +834,11 @@ contains
! We need distribcell if any tallies have distribcell filters.
do i = 1, n_tallies
do j = 1, tallies(i) % n_filters
if (tallies(i) % filters(j) % type == FILTER_DISTRIBCELL) then
do j = 1, size(tallies(i) % filters)
select type(filt => tallies(i) % filters(j) % obj)
type is (DistribcellFilter)
distribcell_active = .true.
if (size(tallies(i) % 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 select
end do
end do
@ -924,13 +861,12 @@ contains
! Set the number of bins in all distribcell filters.
do i = 1, n_tallies
do j = 1, tallies(i) % n_filters
associate (filt => tallies(i) % filters(j))
if (filt % type == FILTER_DISTRIBCELL) then
! Set the number of bins to the number of instances of the cell.
filt % n_bins = cells(filt % int_bins(1)) % instances
end if
end associate
do j = 1, size(tallies(i) % filters)
select type(filt => tallies(i) % filters(j) % obj)
type is (DistribcellFilter)
! Set the number of bins to the number of instances of the cell.
filt % n_bins = cells(filt % cell) % instances
end select
end do
end do
@ -990,10 +926,11 @@ contains
! List all cells referenced in distribcell filters.
do i = 1, n_tallies
do j = 1, tallies(i) % n_filters
if (tallies(i) % filters(j) % type == FILTER_DISTRIBCELL) then
call cell_list % add(tallies(i) % filters(j) % int_bins(1))
end if
do j = 1, size(tallies(i) % filters)
select type(filt => tallies(i) % filters(j) % obj)
type is (DistribcellFilter)
call cell_list % add(filt % cell)
end select
end do
end do

546
src/input_xml.F90
View file

@ -24,7 +24,8 @@ module input_xml
use stl_vector, only: VectorInt, VectorReal, VectorChar
use string, only: to_lower, to_str, str_to_int, str_to_real, &
starts_with, ends_with, tokenize, split_string
use tally_header, only: TallyObject, TallyFilter
use tally_header, only: TallyObject
use tally_filter
use tally_initialize, only: add_tallies
use xml_interface
@ -2667,7 +2668,7 @@ contains
type(ElemKeyValueCI), pointer :: pair_list
type(TallyObject), pointer :: t
type(RegularMesh), pointer :: m
type(TallyFilter), allocatable :: filters(:) ! temporary filters
type(TallyFilterContainer), allocatable :: filters(:) ! temporary filters
type(Node), pointer :: doc => null()
type(Node), pointer :: node_mesh => null()
type(Node), pointer :: node_tal => null()
@ -2907,115 +2908,122 @@ contains
call get_node_list(node_tal, "filter", node_filt_list)
n_filters = get_list_size(node_filt_list)
if (n_filters /= 0) then
! Allocate filters array
allocate(t % filters(n_filters))
! Allocate filters array
t % n_filters = n_filters
allocate(t % filters(n_filters))
READ_FILTERS: do j = 1, n_filters
! Get pointer to filter xml node
call get_list_item(node_filt_list, j, node_filt)
READ_FILTERS: do j = 1, n_filters
! Get pointer to filter xml node
call get_list_item(node_filt_list, j, node_filt)
! Convert filter type to lower case
temp_str = ''
if (check_for_node(node_filt, "type")) &
call get_node_value(node_filt, "type", temp_str)
temp_str = to_lower(temp_str)
! Convert filter type to lower case
temp_str = ''
if (check_for_node(node_filt, "type")) &
call get_node_value(node_filt, "type", temp_str)
temp_str = to_lower(temp_str)
! Determine number of bins
if (check_for_node(node_filt, "bins")) then
if (temp_str == 'energy' .or. temp_str == 'energyout' .or. &
temp_str == 'mu' .or. temp_str == 'polar' .or. &
temp_str == 'azimuthal') then
n_words = get_arraysize_double(node_filt, "bins")
else
n_words = get_arraysize_integer(node_filt, "bins")
end if
! Determine number of bins
if (check_for_node(node_filt, "bins")) then
if (temp_str == 'energy' .or. temp_str == 'energyout' .or. &
temp_str == 'mu' .or. temp_str == 'polar' .or. &
temp_str == 'azimuthal') then
n_words = get_arraysize_double(node_filt, "bins")
else
call fatal_error("Bins not set in filter on tally " &
// trim(to_str(t % id)))
n_words = get_arraysize_integer(node_filt, "bins")
end if
else
call fatal_error("Bins not set in filter on tally " &
// trim(to_str(t % id)))
end if
! Determine type of filter
select case (temp_str)
! 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
case ('distribcell')
! Allocate and declare the filter type
allocate(DistribcellFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (DistribcellFilter)
if (n_words /= 1) call fatal_error("Only one cell can be &
&specified per distribcell filter.")
! Store bins
call get_node_value(node_filt, "bins", filt % cell)
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_DISTRIBCELL) = j
case ('cell')
! Allocate and declare the filter type
allocate(CellFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (CellFilter)
! 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
! Set number of bins
t % filters(j) % n_bins = n_words
filt % n_bins = n_words
allocate(filt % cells(n_words))
call get_node_array(node_filt, "bins", filt % cells)
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_CELL) = j
case ('cellborn')
! Allocate and declare the filter type
allocate(CellbornFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (CellbornFilter)
! 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 ('cellborn')
! Set type of filter
t % filters(j) % type = FILTER_CELLBORN
! Set number of bins
t % filters(j) % n_bins = n_words
filt % n_bins = n_words
allocate(filt % cells(n_words))
call get_node_array(node_filt, "bins", filt % cells)
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_CELLBORN) = j
case ('material')
! Allocate and declare the filter type
allocate(MaterialFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (MaterialFilter)
! 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 ('material')
! Set type of filter
t % filters(j) % type = FILTER_MATERIAL
! Set number of bins
t % filters(j) % n_bins = n_words
filt % n_bins = n_words
allocate(filt % materials(n_words))
call get_node_array(node_filt, "bins", filt % materials)
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_MATERIAL) = j
case ('universe')
! Allocate and declare the filter type
allocate(UniverseFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (UniverseFilter)
! 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 ('universe')
! Set type of filter
t % filters(j) % type = FILTER_UNIVERSE
! Set number of bins
t % filters(j) % n_bins = n_words
filt % n_bins = n_words
allocate(filt % universes(n_words))
call get_node_array(node_filt, "bins", filt % universes)
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_UNIVERSE) = j
case ('surface')
call fatal_error("Surface filter is not yet supported!")
! Allocate and declare the filter type
allocate(SurfaceFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (SurfaceFilter)
! Allocate and store bins
allocate(t % filters(j) % int_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % int_bins)
filt % n_bins = n_words
allocate(filt % surfaces(n_words))
call get_node_array(node_filt, "bins", filt % surfaces)
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_SURFACE) = j
case ('surface')
call fatal_error("Surface filter is not yet supported!")
! Set type of filter
t % filters(j) % type = FILTER_SURFACE
! Set number of bins
t % filters(j) % n_bins = n_words
! 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 ('mesh')
! Set type of filter
t % filters(j) % type = FILTER_MESH
! Check to make sure multiple meshes weren't given
if (n_words /= 1) then
call fatal_error("Can only have one mesh filter specified.")
end if
case ('mesh')
! Allocate and declare the filter type
allocate(MeshFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (MeshFilter)
if (n_words /= 1) call fatal_error("Only one mesh can be &
&specified per mesh filter.")
! Determine id of mesh
call get_node_value(node_filt, "bins", id)
@ -3032,214 +3040,220 @@ contains
! Determine number of bins -- this is assuming that the tally is
! a volume tally and not a surface current tally. If it is a
! surface current tally, the number of bins will get reset later
t % filters(j) % n_bins = product(m % dimension)
filt % n_bins = product(m % dimension)
! Allocate and store index of mesh
allocate(t % filters(j) % int_bins(1))
t % filters(j) % int_bins(1) = i_mesh
case ('energy')
! Set type of filter
t % filters(j) % type = FILTER_ENERGYIN
! Set number of bins
t % filters(j) % n_bins = n_words - 1
! Store the index of the mesh
filt % mesh = i_mesh
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_MESH) = j
case ('energy')
! Allocate and declare the filter type
allocate(EnergyFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (EnergyFilter)
! Allocate and store bins
allocate(t % filters(j) % real_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % real_bins)
filt % n_bins = n_words - 1
allocate(filt % bins(n_words))
call get_node_array(node_filt, "bins", filt % bins)
! We can save tallying time if we know that the tally bins
! match the energy group structure. In that case, the matching bin
! We can save tallying time if we know that the tally bins match
! the energy group structure. In that case, the matching bin
! index is simply the group (after flipping for the different
! ordering of the library and tallying systems).
if (.not. run_CE) then
if (n_words == energy_groups + 1) then
if (all(t % filters(j) % real_bins == &
energy_bins(energy_groups + 1:1:-1))) &
t % energy_matches_groups = .true.
if (all(filt % bins == energy_bins(energy_groups + 1:1:-1))) &
then
filt % matches_transport_groups = .true.
end if
end if
end if
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_ENERGYIN) = j
case ('energyout')
! Set type of filter
t % filters(j) % type = FILTER_ENERGYOUT
! Set number of bins
t % filters(j) % n_bins = n_words - 1
case ('energyout')
! Allocate and declare the filter type
allocate(EnergyoutFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (EnergyoutFilter)
! Allocate and store bins
allocate(t % filters(j) % real_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % real_bins)
filt % n_bins = n_words - 1
allocate(filt % bins(n_words))
call get_node_array(node_filt, "bins", filt % bins)
! We can save tallying time if we know that the tally bins
! match the energy group structure. In that case, the matching bin
! We can save tallying time if we know that the tally bins match
! the energy group structure. In that case, the matching bin
! index is simply the group (after flipping for the different
! ordering of the library and tallying systems).
if (.not. run_CE) then
if (n_words == energy_groups + 1) then
if (all(t % filters(j) % real_bins == &
energy_bins(energy_groups + 1:1:-1))) &
t % energyout_matches_groups = .true.
if (all(filt % bins == energy_bins(energy_groups + 1:1:-1))) &
then
filt % matches_transport_groups = .true.
end if
end if
end if
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_ENERGYOUT) = j
! Set to analog estimator
t % estimator = ESTIMATOR_ANALOG
! Set to analog estimator
t % estimator = ESTIMATOR_ANALOG
case ('delayedgroup')
! Check to see if running in MG mode, because if so, the current
! system isnt set up yet to support delayed group data and thus
! these tallies
if (.not. run_CE) then
call fatal_error("delayedgroup filter on tally " &
// trim(to_str(t % id)) // " not yet supported&
& for multi-group mode.")
end if
! Set type of filter
t % filters(j) % type = FILTER_DELAYEDGROUP
! Set number of bins
t % filters(j) % n_bins = n_words
case ('delayedgroup')
! Check to see if running in MG mode, because if so, the current
! system isnt set up yet to support delayed group data and thus
! these tallies
if (.not. run_CE) then
call fatal_error("delayedgroup filter on tally " &
// trim(to_str(t % id)) // " not yet supported&
& for multi-group mode.")
end if
! Allocate and declare the filter type
allocate(DelayedGroupFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (DelayedGroupFilter)
! Allocate and store bins
allocate(t % filters(j) % int_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % int_bins)
filt % n_bins = n_words
allocate(filt % groups(n_words))
call get_node_array(node_filt, "bins", filt % groups)
! Check bins to make sure all are between 1 and MAX_DELAYED_GROUPS
! Check that bins are all are between 1 and MAX_DELAYED_GROUPS
do d = 1, n_words
if (t % filters(j) % int_bins(d) < 1 .or. &
t % filters(j) % int_bins(d) > MAX_DELAYED_GROUPS) then
if (filt % groups(d) < 1 .or. &
filt % groups(d) > MAX_DELAYED_GROUPS) then
call fatal_error("Encountered delayedgroup bin with index " &
// trim(to_str(t % filters(j) % int_bins(d))) // " that is&
& outside the range of 1 to MAX_DELAYED_GROUPS ( " &
// trim(to_str(filt % groups(d))) // " that is outside &
&the range of 1 to MAX_DELAYED_GROUPS ( " &
// trim(to_str(MAX_DELAYED_GROUPS)) // ")")
end if
end do
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_DELAYEDGROUP) = j
case ('mu')
! Set type of filter
t % filters(j) % type = FILTER_MU
! Set number of bins
t % filters(j) % n_bins = n_words - 1
case ('mu')
! Allocate and declare the filter type
allocate(MuFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (MuFilter)
! Allocate and store bins
allocate(t % filters(j) % real_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % real_bins)
filt % n_bins = n_words - 1
allocate(filt % bins(n_words))
call get_node_array(node_filt, "bins", filt % bins)
! Allow a user to input a lone number which will mean that
! you subivide [-1,1] evenly with the input being the number of bins
! Allow a user to input a lone number which will mean that you
! subdivide [-1,1] evenly with the input being the number of bins
if (n_words == 1) then
Nangle = int(t % filters(j) % real_bins(1))
Nangle = int(filt % bins(1))
if (Nangle > 1) then
t % filters(j) % n_bins = Nangle
filt % n_bins = Nangle
dangle = TWO / real(Nangle,8)
deallocate(t % filters(j) % real_bins)
allocate(t % filters(j) % real_bins(Nangle + 1))
deallocate(filt % bins)
allocate(filt % bins(Nangle + 1))
do iangle = 1, Nangle
t % filters(j) % real_bins(iangle) = -ONE + (iangle - 1) * dangle
filt % bins(iangle) = -ONE + (iangle - 1) * dangle
end do
t % filters(j) % real_bins(Nangle + 1) = ONE
filt % bins(Nangle + 1) = ONE
else
call fatal_error("Number of bins for mu filter must be&
& greater than 1 on tally " // trim(to_str(t % id)) // ".")
& greater than 1 on tally " &
// trim(to_str(t % id)) // ".")
end if
end if
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_MU) = j
! Set to analog estimator
t % estimator = ESTIMATOR_ANALOG
case ('polar')
! Set type of filter
t % filters(j) % type = FILTER_POLAR
! Set number of bins
t % filters(j) % n_bins = n_words - 1
! Set to analog estimator
t % estimator = ESTIMATOR_ANALOG
case ('polar')
! Allocate and declare the filter type
allocate(PolarFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (PolarFilter)
! Allocate and store bins
allocate(t % filters(j) % real_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % real_bins)
filt % n_bins = n_words - 1
allocate(filt % bins(n_words))
call get_node_array(node_filt, "bins", filt % bins)
! Allow a user to input a lone number which will mean that
! you subivide [0,pi] evenly with the input being the number of bins
! Allow a user to input a lone number which will mean that you
! subdivide [0,pi] evenly with the input being the number of bins
if (n_words == 1) then
Nangle = int(t % filters(j) % real_bins(1))
Nangle = int(filt % bins(1))
if (Nangle > 1) then
t % filters(j) % n_bins = Nangle
filt % n_bins = Nangle
dangle = PI / real(Nangle,8)
deallocate(t % filters(j) % real_bins)
allocate(t % filters(j) % real_bins(Nangle + 1))
deallocate(filt % bins)
allocate(filt % bins(Nangle + 1))
do iangle = 1, Nangle
t % filters(j) % real_bins(iangle) = (iangle - 1) * dangle
filt % bins(iangle) = (iangle - 1) * dangle
end do
t % filters(j) % real_bins(Nangle + 1) = PI
filt % bins(Nangle + 1) = PI
else
call fatal_error("Number of bins for polar filter must be&
& greater than 1 on tally " // trim(to_str(t % id)) // ".")
call fatal_error("Number of bins for mu filter must be&
& greater than 1 on tally " &
// trim(to_str(t % id)) // ".")
end if
end if
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_POLAR) = j
case ('azimuthal')
! Set type of filter
t % filters(j) % type = FILTER_AZIMUTHAL
! Set number of bins
t % filters(j) % n_bins = n_words - 1
case ('azimuthal')
! Allocate and declare the filter type
allocate(AzimuthalFilter::t % filters(j) % obj)
select type (filt => t % filters(j) % obj)
type is (AzimuthalFilter)
! Allocate and store bins
allocate(t % filters(j) % real_bins(n_words))
call get_node_array(node_filt, "bins", t % filters(j) % real_bins)
filt % n_bins = n_words - 1
allocate(filt % bins(n_words))
call get_node_array(node_filt, "bins", filt % bins)
! Allow a user to input a lone number which will mean that
! you sub-divide [-pi,pi) evenly with the input being the number of
! Allow a user to input a lone number which will mean that you
! subdivide [-pi,pi) evenly with the input being the number of
! bins
if (n_words == 1) then
Nangle = int(t % filters(j) % real_bins(1))
Nangle = int(filt % bins(1))
if (Nangle > 1) then
t % filters(j) % n_bins = Nangle
filt % n_bins = Nangle
dangle = TWO * PI / real(Nangle,8)
deallocate(t % filters(j) % real_bins)
allocate(t % filters(j) % real_bins(Nangle + 1))
deallocate(filt % bins)
allocate(filt % bins(Nangle + 1))
do iangle = 1, Nangle
t % filters(j) % real_bins(iangle) = -PI + (iangle - 1) * dangle
filt % bins(iangle) = -PI + (iangle - 1) * dangle
end do
t % filters(j) % real_bins(Nangle + 1) = PI
filt % bins(Nangle + 1) = PI
else
call fatal_error("Number of bins for azimuthal filter must be&
& greater than 1 on tally " // trim(to_str(t % id)) // ".")
call fatal_error("Number of bins for mu filter must be&
& greater than 1 on tally " &
// trim(to_str(t % id)) // ".")
end if
end if
case default
! Specified tally filter is invalid, raise error
call fatal_error("Unknown filter type '" &
// trim(temp_str) // "' on tally " &
// trim(to_str(t % id)) // ".")
end select
! Set the filter index in the tally find_filter array
t % find_filter(FILTER_AZIMUTHAL) = j
! Set find_filter, e.g. if filter(3) has type FILTER_CELL, then
! find_filter(FILTER_CELL) would be set to 3.
case default
! Specified tally filter is invalid, raise error
call fatal_error("Unknown filter type '" &
// trim(temp_str) // "' on tally " &
// trim(to_str(t % id)) // ".")
t % find_filter(t % filters(j) % type) = j
end select
end do READ_FILTERS
end do READ_FILTERS
! Check that both cell and surface weren't specified
if (t % find_filter(FILTER_CELL) > 0 .and. &
t % find_filter(FILTER_SURFACE) > 0) then
call fatal_error("Cannot specify both cell and surface filters for &
&tally " // trim(to_str(t % id)))
end if
else
! No filters were specified
t % n_filters = 0
! Check that both cell and surface weren't specified
if (t % find_filter(FILTER_CELL) > 0 .and. &
t % find_filter(FILTER_SURFACE) > 0) then
call fatal_error("Cannot specify both cell and surface filters for &
&tally " // trim(to_str(t % id)))
end if
! =======================================================================
@ -3275,8 +3289,9 @@ contains
if (trim(sarray(j)) == 'total') then
! Check if a delayedgroup filter is present for this tally
do l = 1, t % n_filters
if (t % filters(l) % type == FILTER_DELAYEDGROUP) then
do l = 1, size(t % filters)
select type(filt => t % filters(l) % obj)
type is (DelayedGroupFilter)
call warning("A delayedgroup filter was used on a total &
&nuclide tally. Cross section libraries are not &
&guaranteed to have the same delayed group structure &
@ -3285,7 +3300,7 @@ contains
&all isotopes while the JEFF 3.1.1 library has the same &
&delayed group structure across all isotopes. Use with &
&caution!")
end if
end select
end do
t % nuclide_bins(j) = -1
@ -3334,8 +3349,9 @@ contains
t % n_nuclide_bins = 1
! Check if a delayedgroup filter is present for this tally
do l = 1, t % n_filters
if (t % filters(l) % type == FILTER_DELAYEDGROUP) then
do l = 1, size(t % filters)
select type(filt => t % filters(l) % obj)
type is (DelayedGroupFilter)
call warning("A delayedgroup filter was used on a total nuclide &
&tally. Cross section libraries are not guaranteed to have the&
& same delayed group structure across all isotopes. In &
@ -3343,7 +3359,7 @@ contains
&group structure across all isotopes while the JEFF 3.1.1 &
&library has the same delayed group structure across all &
&isotopes. Use with caution!")
end if
end select
end do
end if
@ -3659,36 +3675,42 @@ contains
&filter.")
end if
! Get pointer to mesh
i_mesh = t % filters(k) % int_bins(1)
m => meshes(i_mesh)
! Declare the type of the mesh filter
select type(filt => t % filters(k) % obj)
type is (MeshFilter)
! We need to increase the dimension by one since we also need
! currents coming into and out of the boundary mesh cells.
t % filters(k) % n_bins = product(m % dimension + 1)
! Get pointer to mesh
i_mesh = filt % mesh
m => meshes(i_mesh)
! We need to increase the dimension by one since we also need
! currents coming into and out of the boundary mesh cells.
filt % n_bins = product(m % dimension + 1)
end select
! Copy filters to temporary array
allocate(filters(t % n_filters + 1))
filters(1:t % n_filters) = t % filters
allocate(filters(size(t % filters) + 1))
filters(1:size(t % filters)) = t % filters
! Move allocation back -- filters becomes deallocated during
! this call
call move_alloc(FROM=filters, TO=t%filters)
! Add surface filter
t % n_filters = t % n_filters + 1
t % filters(t % n_filters) % type = FILTER_SURFACE
t % filters(t % n_filters) % n_bins = 2 * m % n_dimension
allocate(t % filters(t % n_filters) % int_bins(&
2 * m % n_dimension))
if (m % n_dimension == 2) then
t % filters(t % n_filters) % int_bins = (/ IN_RIGHT, &
OUT_RIGHT, IN_FRONT, OUT_FRONT /)
elseif (m % n_dimension == 3) then
t % filters(t % n_filters) % int_bins = (/ IN_RIGHT, &
OUT_RIGHT, IN_FRONT, OUT_FRONT, IN_TOP, OUT_TOP /)
end if
t % find_filter(FILTER_SURFACE) = t % n_filters
n_filters = size(t % filters)
allocate(SurfaceFilter :: t % filters(n_filters) % obj)
select type (filt => t % filters(size(t % filters)) % obj)
type is (SurfaceFilter)
filt % n_bins = 2 * m % n_dimension
allocate(filt % surfaces(2 * m % n_dimension))
if (m % n_dimension == 2) then
filt % surfaces = (/ IN_RIGHT, OUT_RIGHT, IN_FRONT, OUT_FRONT /)
elseif (m % n_dimension == 3) then
filt % surfaces = (/ IN_RIGHT, OUT_RIGHT, IN_FRONT, OUT_FRONT,&
IN_TOP, OUT_TOP /)
end if
end select
t % find_filter(FILTER_SURFACE) = size(t % filters)
case ('events')
t % score_bins(j) = SCORE_EVENTS
@ -4398,9 +4420,11 @@ contains
&meshlines on plot " // trim(to_str(pl % id)))
end if
i_mesh = cmfd_tallies(1) % &
filters(cmfd_tallies(1) % find_filter(FILTER_MESH)) % &
int_bins(1)
select type(filt => cmfd_tallies(1) % &
filters(cmfd_tallies(1) % find_filter(FILTER_MESH)) % obj)
type is (MeshFilter)
i_mesh = filt % mesh
end select
pl % meshlines_mesh => meshes(i_mesh)
case ('entropy')

112
src/mesh.F90
View file

@ -290,34 +290,42 @@ contains
! Check if line intersects left surface -- calculate the intersection point
! y
yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then
yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
end if
! Check if line intersects back surface -- calculate the intersection point
! x
xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1) then
intersects = .true.
return
if ((y0 < ym0 .and. y1 > ym0) .or. (y0 > ym0 .and. y1 < ym0)) then
xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1) then
intersects = .true.
return
end if
end if
! Check if line intersects right surface -- calculate the intersection
! point y
yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then
yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
end if
! Check if line intersects front surface -- calculate the intersection point
! x
xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1) then
intersects = .true.
return
if ((y0 < ym1 .and. y1 > ym1) .or. (y0 > ym1 .and. y1 < ym1)) then
xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1) then
intersects = .true.
return
end if
end if
end function mesh_intersects_2d
@ -359,56 +367,68 @@ contains
! Check if line intersects left surface -- calculate the intersection point
! (y,z)
yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then
yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
end if
! Check if line intersects back surface -- calculate the intersection point
! (x,z)
xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
if ((y0 < ym0 .and. y1 > ym0) .or. (y0 > ym0 .and. y1 < ym0)) then
xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
end if
! Check if line intersects bottom surface -- calculate the intersection
! point (x,y)
xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0)
if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
if ((z0 < zm0 .and. z1 > zm0) .or. (z0 > zm0 .and. z1 < zm0)) then
xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0)
if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
end if
! Check if line intersects right surface -- calculate the intersection point
! (y,z)
yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then
yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0)
zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0)
if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
end if
! Check if line intersects front surface -- calculate the intersection point
! (x,z)
xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
if ((y0 < ym1 .and. y1 > ym1) .or. (y0 > ym1 .and. y1 < ym1)) then
xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0)
zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0)
if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
intersects = .true.
return
end if
end if
! Check if line intersects top surface -- calculate the intersection point
! (x,y)
xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0)
if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
if ((z0 < zm1 .and. z1 > zm1) .or. (z0 > zm1 .and. z1 < zm1)) then
xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0)
yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0)
if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
intersects = .true.
return
end if
end if
end function mesh_intersects_3d

422
src/output.F90
View file

@ -17,6 +17,7 @@ module output
use sab_header, only: SAlphaBeta
use string, only: to_upper, to_str
use tally_header, only: TallyObject
use tally_filter
implicit none
@ -736,7 +737,6 @@ contains
integer :: k ! loop index for scoring bins
integer :: n ! loop index for nuclides
integer :: l ! loop index for user scores
integer :: type ! type of tally filter
integer :: indent ! number of spaces to preceed output
integer :: filter_index ! index in results array for filters
integer :: score_index ! scoring bin index
@ -747,7 +747,6 @@ 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(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
@ -756,21 +755,6 @@ contains
! Skip if there are no tallies
if (n_tallies == 0) return
! Initialize names for tally filter types
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"
filter_name(FILTER_MU) = "Change-in-Angle"
filter_name(FILTER_POLAR) = "Polar Angle"
filter_name(FILTER_AZIMUTHAL) = "Azimuthal Angle"
filter_name(FILTER_DELAYEDGROUP) = "Delayed Group"
! Initialize names for scores
score_names(abs(SCORE_FLUX)) = "Flux"
score_names(abs(SCORE_TOTAL)) = "Total Reaction Rate"
@ -841,14 +825,14 @@ contains
! to be used for a given tally.
! Initialize bins, filter level, and indentation
matching_bins(1:t%n_filters) = 0
matching_bins(1:size(t % filters)) = 0
j = 1
indent = 0
print_bin: do
find_bin: do
! Check for no filters
if (t % n_filters == 0) exit find_bin
if (size(t % filters) == 0) exit find_bin
! Increment bin combination
matching_bins(j) = matching_bins(j) + 1
@ -856,7 +840,7 @@ contains
! =================================================================
! REACHED END OF BINS FOR THIS FILTER, MOVE TO NEXT FILTER
if (matching_bins(j) > t % filters(j) % n_bins) then
if (matching_bins(j) > t % filters(j) % obj % n_bins) then
! If this is the first filter, then exit
if (j == 1) exit print_bin
@ -869,12 +853,11 @@ contains
else
! Check if this is last filter
if (j == t % n_filters) exit find_bin
if (j == size(t % filters)) exit find_bin
! Print current filter information
type = t % filters(j) % type
write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(filter_name(type)), trim(get_label(t, j))
write(UNIT=unit_tally, FMT='(1X,2A)') repeat(" ", indent), &
trim(t % filters(j) % obj % text_label(matching_bins(j)))
indent = indent + 2
j = j + 1
end if
@ -882,25 +865,25 @@ contains
end do find_bin
! Print filter information
if (t % n_filters > 0) then
type = t % filters(j) % type
write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(filter_name(type)), trim(get_label(t, j))
if (size(t % filters) > 0) then
write(UNIT=unit_tally, FMT='(1X,2A)') repeat(" ", indent), &
trim(t % filters(j) % obj % text_label(matching_bins(j)))
end if
! Determine scoring index for this bin combination -- note that unlike
! in the score_tally subroutine, we have to use max(bins,1) since all
! bins below the lowest filter level will be zeros
if (t % n_filters > 0) then
filter_index = sum((max(matching_bins(1:t%n_filters),1) - 1) * t % stride) + 1
if (size(t % filters) > 0) then
filter_index = sum((max(matching_bins(1:size(t % filters)),1) - 1) &
* t % stride) + 1
else
filter_index = 1
end if
! Write results for this filter bin combination
score_index = 0
if (t % n_filters > 0) indent = indent + 2
if (size(t % filters) > 0) indent = indent + 2
do n = 1, t % n_nuclide_bins
! Write label for nuclide
i_nuclide = t % nuclide_bins(n)
@ -972,7 +955,7 @@ contains
end do
indent = indent - 2
if (t % n_filters == 0) exit print_bin
if (size(t % filters) == 0) exit print_bin
end do print_bin
@ -1009,16 +992,19 @@ contains
! Get pointer to mesh
i_filter_mesh = t % find_filter(FILTER_MESH)
i_filter_surf = t % find_filter(FILTER_SURFACE)
m => meshes(t % filters(i_filter_mesh) % int_bins(1))
select type(filt => t % filters(i_filter_mesh) % obj)
type is (MeshFilter)
m => meshes(filt % mesh)
end select
! initialize bins array
matching_bins(1:t%n_filters) = 1
matching_bins(1:size(t % filters)) = 1
! determine how many energy in bins there are
i_filter_ein = t % find_filter(FILTER_ENERGYIN)
if (i_filter_ein > 0) then
print_ebin = .true.
n = t % filters(i_filter_ein) % n_bins
n = t % filters(i_filter_ein) % obj % n_bins
else
print_ebin = .false.
n = 1
@ -1041,22 +1027,25 @@ contains
matching_bins(i_filter_ein) = l
! Write incoming energy bin
write(UNIT=unit_tally, FMT='(3X,A,1X,A)') &
"Incoming Energy", trim(get_label(t, i_filter_ein))
write(UNIT=unit_tally, FMT='(3X,A)') &
trim(t % filters(i_filter_ein) % obj % text_label( &
matching_bins(i_filter_ein)))
end if
! Left Surface
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Left", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Left", &
to_str(t % results(1,filter_index) % sum), &
@ -1066,14 +1055,16 @@ contains
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Right", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Right", &
to_str(t % results(1,filter_index) % sum), &
@ -1083,14 +1074,16 @@ contains
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Back", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
matching_bins(i_filter_surf) = OUT_FRONT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Back", &
to_str(t % results(1,filter_index) % sum), &
@ -1100,14 +1093,16 @@ contains
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Front", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
matching_bins(i_filter_surf) = OUT_FRONT
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Front", &
to_str(t % results(1,filter_index) % sum), &
@ -1117,14 +1112,16 @@ contains
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.)
matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Bottom", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
matching_bins(i_filter_surf) = OUT_TOP
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Bottom", &
to_str(t % results(1,filter_index) % sum), &
@ -1134,14 +1131,16 @@ contains
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Top", &
to_str(t % results(1,filter_index) % sum), &
trim(to_str(t % results(1,filter_index) % sum_sq))
matching_bins(i_filter_surf) = OUT_TOP
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
filter_index = sum((matching_bins(1:size(t % filters)) - 1) &
* t % stride) + 1
write(UNIT=unit_tally, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Top", &
to_str(t % results(1,filter_index) % sum), &
@ -1154,329 +1153,4 @@ contains
end subroutine write_surface_current
!===============================================================================
! GET_LABEL returns a label for a cell/surface/etc given a tally, filter type,
! and corresponding bin
!===============================================================================
function get_label(t, i_filter) result(label)
type(TallyObject), intent(in) :: t ! tally object
integer, intent(in) :: i_filter ! index in filters array
character(MAX_LINE_LEN) :: 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(RegularMesh), pointer :: m
type(Universe), pointer :: univ
bin = matching_bins(i_filter)
select case(t % filters(i_filter) % type)
case (FILTER_UNIVERSE)
i = t % filters(i_filter) % int_bins(bin)
label = to_str(universes(i) % id)
case (FILTER_MATERIAL)
i = t % filters(i_filter) % int_bins(bin)
label = to_str(materials(i) % id)
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) % int_bins(1), &
univ, bin-1, offset, label)
case (FILTER_SURFACE)
i = t % filters(i_filter) % int_bins(bin)
label = to_str(surfaces(i)%obj%id)
case (FILTER_MESH)
m => meshes(t % filters(i_filter) % int_bins(1))
allocate(ijk(m % n_dimension))
call bin_to_mesh_indices(m, bin, ijk)
if (m % n_dimension == 2) then
label = "Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ")"
elseif (m % n_dimension == 3) then
label = "Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
end if
case (FILTER_ENERGYIN, FILTER_ENERGYOUT, FILTER_MU, FILTER_POLAR, &
FILTER_AZIMUTHAL)
E0 = t % filters(i_filter) % real_bins(bin)
E1 = t % filters(i_filter) % real_bins(bin + 1)
label = "[" // trim(to_str(E0)) // ", " // trim(to_str(E1)) // ")"
case (FILTER_DELAYEDGROUP)
i = t % filters(i_filter) % int_bins(bin)
label = to_str(i)
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(goal, univ, final, offset, path)
integer, intent(in) :: goal ! The target cell index
type(Universe), intent(in) :: univ ! Universe to begin search
integer, intent(in) :: final ! Target offset
integer, intent(inout) :: offset ! Current offset
character(*), intent(inout) :: path ! Path to offset
integer :: map ! Index in maps vector
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
! Get the distribcell index for this cell
map = cells(goal) % distribcell_index
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
! If the cell matches the goal and the offset matches final, write to the
! geometry stack
if (univ % cells(i) == goal .and. offset == final) then
c => cells(univ % cells(i))
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(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(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(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
if (k == lat % n_rings .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 - lat % n_rings)) // "," // &
trim(to_str(l - lat % n_rings)) // "," // &
trim(to_str(m)) // ")"
call find_offset(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 - lat % n_rings)) // "," // &
trim(to_str(old_l - lat % n_rings)) // "," // &
trim(to_str(old_m)) // ")"
call find_offset(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

34
src/particle_header.F90
View file

@ -59,10 +59,13 @@ module particle_header
logical :: alive ! is particle alive?
! Pre-collision physical data
real(8) :: last_xyz(3) ! previous coordinates
real(8) :: last_uvw(3) ! previous direction coordinates
real(8) :: last_wgt ! pre-collision particle weight
real(8) :: absorb_wgt ! weight absorbed for survival biasing
real(8) :: last_xyz_current(3) ! coordinates of the last collision or
! reflective/periodic surface crossing
! for current tallies
real(8) :: last_xyz(3) ! previous coordinates
real(8) :: last_uvw(3) ! previous direction coordinates
real(8) :: last_wgt ! pre-collision particle weight
real(8) :: absorb_wgt ! weight absorbed for survival biasing
! What event last took place
logical :: fission ! did the particle cause implicit fission
@ -193,20 +196,21 @@ contains
call this % initialize()
! copy attributes from source bank site
this % wgt = src % wgt
this % last_wgt = src % wgt
this % coord(1) % xyz = src % xyz
this % coord(1) % uvw = src % uvw
this % last_xyz = src % xyz
this % last_uvw = src % uvw
this % wgt = src % wgt
this % last_wgt = src % wgt
this % coord(1) % xyz = src % xyz
this % coord(1) % uvw = src % uvw
this % last_xyz_current = src % xyz
this % last_xyz = src % xyz
this % last_uvw = src % uvw
if (run_CE) then
this % E = src % E
this % E = src % E
else
this % g = int(src % E)
this % last_g = int(src % E)
this % E = energy_bin_avg(this % g)
this % g = int(src % E)
this % last_g = int(src % E)
this % E = energy_bin_avg(this % g)
end if
this % last_E = this % E
this % last_E = this % E
end subroutine initialize_from_source

11
src/particle_restart.F90
View file

@ -107,11 +107,12 @@ contains
end if
! Set particle last attributes
p % last_wgt = p % wgt
p % last_xyz = p % coord(1)%xyz
p % last_uvw = p % coord(1)%uvw
p % last_E = p % E
p % last_g = p % g
p % last_wgt = p % wgt
p % last_xyz_current = p % coord(1)%xyz
p % last_xyz = p % coord(1)%xyz
p % last_uvw = p % coord(1)%uvw
p % last_E = p % E
p % last_g = p % g
! Close hdf5 file
call file_close(file_id)

50
src/state_point.F90
View file

@ -237,57 +237,13 @@ contains
end select
call write_dataset(tally_group, "n_realizations", &
tally % n_realizations)
call write_dataset(tally_group, "n_filters", tally % n_filters)
call write_dataset(tally_group, "n_filters", size(tally % filters))
! Write filter information
FILTER_LOOP: do j = 1, tally % n_filters
FILTER_LOOP: do j = 1, size(tally % filters)
filter_group = create_group(tally_group, "filter " // &
trim(to_str(j)))
! Write name of type
select case (tally % filters(j) % type)
case(FILTER_UNIVERSE)
call write_dataset(filter_group, "type", "universe")
case(FILTER_MATERIAL)
call write_dataset(filter_group, "type", "material")
case(FILTER_CELL)
call write_dataset(filter_group, "type", "cell")
case(FILTER_CELLBORN)
call write_dataset(filter_group, "type", "cellborn")
case(FILTER_SURFACE)
call write_dataset(filter_group, "type", "surface")
case(FILTER_MESH)
call write_dataset(filter_group, "type", "mesh")
case(FILTER_ENERGYIN)
call write_dataset(filter_group, "type", "energy")
case(FILTER_ENERGYOUT)
call write_dataset(filter_group, "type", "energyout")
case(FILTER_MU)
call write_dataset(filter_group, "type", "mu")
case(FILTER_POLAR)
call write_dataset(filter_group, "type", "polar")
case(FILTER_AZIMUTHAL)
call write_dataset(filter_group, "type", "azimuthal")
case(FILTER_DISTRIBCELL)
call write_dataset(filter_group, "type", "distribcell")
case(FILTER_DELAYEDGROUP)
call write_dataset(filter_group, "type", "delayedgroup")
end select
call write_dataset(filter_group, "n_bins", &
tally % filters(j) % n_bins)
if (tally % filters(j) % type == FILTER_ENERGYIN .or. &
tally % filters(j) % type == FILTER_ENERGYOUT .or. &
tally % filters(j) % type == FILTER_MU .or. &
tally % filters(j) % type == FILTER_POLAR .or. &
tally % filters(j) % type == FILTER_AZIMUTHAL) then
call write_dataset(filter_group, "bins", &
tally % filters(j) % real_bins)
else
call write_dataset(filter_group, "bins", &
tally % filters(j) % int_bins)
end if
call tally % filters(j) % obj % to_statepoint(filter_group)
call close_group(filter_group)
end do FILTER_LOOP

76
src/summary.F90
View file

@ -3,7 +3,7 @@ module summary
use constants
use endf, only: reaction_name
use geometry_header, only: Cell, Universe, Lattice, RectLattice, &
&HexLattice, BASE_UNIVERSE
&HexLattice
use global
use hdf5_interface
use material_header, only: Material
@ -13,7 +13,6 @@ module summary
use surface_header
use string, only: to_str
use tally_header, only: TallyObject
use output, only: find_offset
use hdf5
@ -584,10 +583,6 @@ contains
type(RegularMesh), pointer :: m
type(TallyObject), pointer :: t
integer :: offset ! distibcell offset
character(MAX_LINE_LEN), allocatable :: paths(:) ! distribcell paths array
character(MAX_LINE_LEN) :: path ! distribcell path
tallies_group = create_group(file_id, "tallies")
! Write total number of meshes
@ -628,74 +623,11 @@ contains
call write_dataset(tally_group, "name", t%name)
! Write number of filters
call write_dataset(tally_group, "n_filters", t%n_filters)
call write_dataset(tally_group, "n_filters", size(t % filters))
FILTER_LOOP: do j = 1, t % n_filters
FILTER_LOOP: do j = 1, size(t % filters)
filter_group = create_group(tally_group, "filter " // trim(to_str(j)))
! Write number of bins for this filter
call write_dataset(filter_group, "n_bins", t % filters(j) % n_bins)
! Write filter bins
if (t % filters(j) % type == FILTER_ENERGYIN .or. &
t % filters(j)% type == FILTER_ENERGYOUT .or. &
t % filters(j) % type == FILTER_MU .or. &
t % filters(j) % type == FILTER_POLAR .or. &
t % filters(j) % type == FILTER_AZIMUTHAL) then
call write_dataset(filter_group, "bins", t % filters(j) % real_bins)
else
call write_dataset(filter_group, "bins", t % filters(j) % int_bins)
end if
! Write paths to reach each distribcell instance
if (t % filters(j) % type == FILTER_DISTRIBCELL) then
! Allocate array of strings for each distribcell path
allocate(paths(t % filters(j) % n_bins))
! Store path for each distribcell instance
do k = 1, t % filters(j) % n_bins
path = ''
offset = 1
call find_offset(t % filters(j) % int_bins(1), &
universes(BASE_UNIVERSE), k, offset, path)
paths(k) = path
end do
! Write array of distribcell paths to summary file
call write_dataset(filter_group, "paths", paths)
deallocate(paths)
end if
! Write name of type
select case (t%filters(j)%type)
case(FILTER_UNIVERSE)
call write_dataset(filter_group, "type", "universe")
case(FILTER_MATERIAL)
call write_dataset(filter_group, "type", "material")
case(FILTER_CELL)
call write_dataset(filter_group, "type", "cell")
case(FILTER_CELLBORN)
call write_dataset(filter_group, "type", "cellborn")
case(FILTER_SURFACE)
call write_dataset(filter_group, "type", "surface")
case(FILTER_MESH)
call write_dataset(filter_group, "type", "mesh")
case(FILTER_ENERGYIN)
call write_dataset(filter_group, "type", "energy")
case(FILTER_ENERGYOUT)
call write_dataset(filter_group, "type", "energyout")
case(FILTER_DISTRIBCELL)
call write_dataset(filter_group, "type", "distribcell")
case(FILTER_MU)
call write_dataset(filter_group, "type", "mu")
case(FILTER_POLAR)
call write_dataset(filter_group, "type", "polar")
case(FILTER_AZIMUTHAL)
call write_dataset(filter_group, "type", "azimuthal")
case(FILTER_DELAYEDGROUP)
call write_dataset(filter_group, "type", "delayedgroup")
end select
call t % filters(j) % obj % to_summary(filter_group)
call close_group(filter_group)
end do FILTER_LOOP

1806
src/tally.F90
View file

File diff suppressed because it is too large Load diff

1491
src/tally_filter.F90 Normal file
View file

File diff suppressed because it is too large Load diff

105
src/tally_filter_header.F90 Normal file
View file

@ -0,0 +1,105 @@
module tally_filter_header
use constants, only: MAX_LINE_LEN
use particle_header, only: Particle
use hdf5
implicit none
!===============================================================================
! TALLYFILTER describes a filter that limits what events score to a tally. For
! example, a cell filter indicates that only particles in a specified cell
! should score to the tally.
!===============================================================================
type, abstract :: TallyFilter
integer :: n_bins = 0
contains
procedure(get_next_bin_), deferred :: get_next_bin
procedure(to_statepoint_), deferred :: to_statepoint
procedure :: to_summary => filter_to_summary
procedure(text_label_), deferred :: text_label
procedure :: initialize => filter_initialize
end type TallyFilter
abstract interface
!===============================================================================
! GET_NEXT_BIN gives the index for the next valid filter bin and a weight that
! will be applied to the flux.
!
! In principle, a filter can have multiple valid bins. If current_bin =
! NO_BIN_FOUND, then this method should give the first valid bin. Providing the
! first valid bin should then give the second valid bin, and so on. When there
! are no valid bins left, the next_bin should be NO_VALID_BIN.
subroutine get_next_bin_(this, p, estimator, current_bin, next_bin, weight)
import TallyFilter
import Particle
class(TallyFilter), intent(in) :: this
type(Particle), intent(in) :: p
integer, intent(in) :: estimator
integer, value, intent(in) :: current_bin
integer, intent(out) :: next_bin
real(8), intent(out) :: weight
end subroutine get_next_bin_
!===============================================================================
! TO_STATEPOINT writes all the information needed to reconstruct the filter to
! the given filter_group.
subroutine to_statepoint_(this, filter_group)
import TallyFilter
import HID_T
class(TallyFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
end subroutine to_statepoint_
!===============================================================================
! TEXT_LABEL returns a string describing the given filter bin. For example, an
! energy filter might return the string "Incoming Energy [0.625E-6, 20.0)".
! This is used to write the tallies.out file.
function text_label_(this, bin) result(label)
import TallyFilter
import MAX_LINE_LEN
class(TallyFilter), intent(in) :: this
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
end function text_label_
end interface
!===============================================================================
! TALLYFILTERCONTAINER contains an allocatable TallyFilter object for arrays of
! TallyFilters
!===============================================================================
type TallyFilterContainer
class(TallyFilter), allocatable :: obj
end type TallyFilterContainer
contains
!===============================================================================
! TO_SUMMARY writes all the information needed to reconstruct the filter to the
! given filter_group. If this procedure is not overridden by the derived class,
! then it will call to_statepoint by default.
subroutine filter_to_summary(this, filter_group)
class(TallyFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
call this % to_statepoint(filter_group)
end subroutine filter_to_summary
!===============================================================================
! INITIALIZE sets up any internal data, as necessary. If this procedure is not
! overriden by the derived class, then it will do nothing by default.
subroutine filter_initialize(this)
class(TallyFilter), intent(inout) :: this
end subroutine filter_initialize
end module tally_filter_header

59
src/tally_header.F90
View file

@ -1,41 +1,13 @@
module tally_header
use constants, only: NONE, N_FILTER_TYPES
use trigger_header, only: TriggerObject
use constants, only: NONE, N_FILTER_TYPES
use tally_filter_header, only: TallyFilterContainer
use trigger_header, only: TriggerObject
use, intrinsic :: ISO_C_BINDING
implicit none
!===============================================================================
! TALLYMAPELEMENT gives an index to a tally which is to be scored and the
! corresponding bin for the filter variable
!===============================================================================
type TallyMapElement
integer :: index_tally
integer :: index_bin
end type TallyMapElement
!===============================================================================
! TALLYMAPITEM contains a list of tally/bin combinations for each mappable
! filter bin specified.
!===============================================================================
type TallyMapItem
type(TallyMapElement), allocatable :: elements(:)
end type TallyMapItem
!===============================================================================
! TALLYMAP contains a list of pairs of indices to tallies and the corresponding
! bin for a given filter. There is one TallyMap for each mappable filter
! type. The items array is as long as the corresponding array for that filter,
! e.g. for tally_maps(FILTER_CELL), items is n_cells long.
!===============================================================================
type TallyMap
type(TallyMapItem), allocatable :: items(:)
end type TallyMap
!===============================================================================
! TALLYRESULT provides accumulation of results in a particular tally bin
!===============================================================================
@ -46,19 +18,6 @@ module tally_header
real(C_DOUBLE) :: sum_sq = 0.
end type TallyResult
!===============================================================================
! TALLYFILTER describes a filter that limits what events score to a tally. For
! example, a cell filter indicates that only particles in a specified cell
! should score to the tally.
!===============================================================================
type TallyFilter
integer :: type = NONE
integer :: n_bins = 0
integer, allocatable :: int_bins(:)
real(8), allocatable :: real_bins(:) ! Only used for energy filters
end type TallyFilter
!===============================================================================
! TALLYOBJECT describes a user-specified tally. The region of phase space to
! tally in is given by the TallyFilters and the results are stored in a
@ -73,11 +32,7 @@ module tally_header
integer :: type ! volume, surface current
integer :: estimator ! collision, track-length
real(8) :: volume ! volume of region
! Information about what filters should be used
integer :: n_filters ! Number of filters
type(TallyFilter), allocatable :: filters(:) ! Filter data (type/bins)
type(TallyFilterContainer), allocatable :: filters(:)
! The stride attribute is used for determining the index in the results
! array for a matching_bin combination. Since multiple dimensions are
@ -124,10 +79,6 @@ module tally_header
! Tally precision triggers
integer :: n_triggers = 0 ! # of triggers
type(TriggerObject), allocatable :: triggers(:) ! Array of triggers
! Multi-Group Specific Information To Enable Rapid Tallying
logical :: energy_matches_groups = .false.
logical :: energyout_matches_groups = .false.
end type TallyObject
end module tally_header

103
src/tally_initialize.F90
View file

@ -2,7 +2,7 @@ module tally_initialize
use constants
use global
use tally_header, only: TallyObject, TallyMapElement, TallyMapItem
use tally_header, only: TallyObject
implicit none
private
@ -23,7 +23,6 @@ contains
allocate(global_tallies(N_GLOBAL_TALLIES))
call setup_tally_arrays()
call setup_tally_maps()
end subroutine configure_tallies
@ -45,17 +44,17 @@ contains
t => tallies(i)
! Allocate stride and matching_bins arrays
allocate(t % stride(t % n_filters))
max_n_filters = max(max_n_filters, t % n_filters)
allocate(t % stride(size(t % filters)))
max_n_filters = max(max_n_filters, size(t % filters))
! The filters are traversed in opposite order so that the last filter has
! the shortest stride in memory and the first filter has the largest
! stride
n = 1
STRIDE: do j = t % n_filters, 1, -1
STRIDE: do j = size(t % filters), 1, -1
t % stride(j) = n
n = n * t % filters(j) % n_bins
n = n * t % filters(j) % obj % n_bins
end do STRIDE
! Set total number of filter and scoring bins
@ -70,101 +69,11 @@ contains
! Allocate array for matching filter bins
!$omp parallel
allocate(matching_bins(max_n_filters))
allocate(filter_weights(max_n_filters))
!$omp end parallel
end subroutine setup_tally_arrays
!===============================================================================
! SETUP_TALLY_MAPS creates a map that allows a quick determination of which
! tallies and bins need to be scored to when a particle makes a collision. This
! subroutine also sets the stride attribute for each tally as well as allocating
! storage for the results array.
!===============================================================================
subroutine setup_tally_maps()
integer :: i ! loop index for tallies
integer :: j ! loop index for filters
integer :: k ! loop index for bins
integer :: bin ! filter bin entries
integer :: type ! type of tally filter
type(TallyObject), pointer :: t
! allocate tally map array -- note that we don't need a tally map for the
! energy_in and energy_out filters
allocate(tally_maps(N_FILTER_TYPES - 3))
! allocate list of items for each different filter type
allocate(tally_maps(FILTER_UNIVERSE) % items(n_universes))
allocate(tally_maps(FILTER_MATERIAL) % items(n_materials))
allocate(tally_maps(FILTER_CELL) % items(n_cells))
allocate(tally_maps(FILTER_CELLBORN) % items(n_cells))
allocate(tally_maps(FILTER_SURFACE) % items(n_surfaces))
TALLY_LOOP: do i = 1, n_tallies
! Get pointer to tally
t => tallies(i)
! No need to set up tally maps for surface current tallies
if (t % type == TALLY_SURFACE_CURRENT) cycle
FILTER_LOOP: do j = 1, t % n_filters
! Determine type of filter
type = t % filters(j) % type
if (type == FILTER_CELL .or. type == FILTER_SURFACE .or. &
type == FILTER_MATERIAL .or. type == FILTER_UNIVERSE .or. &
type == FILTER_CELLBORN) then
! Add map elements
BIN_LOOP: do k = 1, t % filters(j) % n_bins
bin = t % filters(j) % int_bins(k)
call add_map_element(tally_maps(type) % items(bin), i, k)
end do BIN_LOOP
end if
end do FILTER_LOOP
end do TALLY_LOOP
end subroutine setup_tally_maps
!===============================================================================
! ADD_MAP_ELEMENT adds a pair of tally and bin indices to the list for a given
! cell/surface/etc.
!===============================================================================
subroutine add_map_element(item, index_tally, index_bin)
type(TallyMapItem), intent(inout) :: item
integer, intent(in) :: index_tally ! index in tallies array
integer, intent(in) :: index_bin ! index in bins array
integer :: n ! size of elements array
type(TallyMapElement), allocatable :: temp(:)
if (.not. allocated(item % elements)) then
allocate(item % elements(1))
item % elements(1) % index_tally = index_tally
item % elements(1) % index_bin = index_bin
else
! determine size of elements array
n = size(item % elements)
! allocate temporary storage and copy elements
allocate(temp(n+1))
temp(1:n) = item % elements
! move allocation back to main array
call move_alloc(FROM=temp, TO=item%elements)
! set new element
item % elements(n+1) % index_tally = index_tally
item % elements(n+1) % index_bin = index_bin
end if
end subroutine add_map_element
!===============================================================================
! ADD_TALLIES extends the tallies array with a new group of tallies and assigns
! pointers to each group. This is called once for user tallies, once for CMFD

5
src/tracking.F90
View file

@ -191,7 +191,7 @@ contains
p % fission = .false.
! Save coordinates for tallying purposes
p % last_xyz = p % coord(1) % xyz
p % last_xyz_current = p % coord(1) % xyz
! Set last material to none since cross sections will need to be
! re-evaluated
@ -211,6 +211,9 @@ contains
end do
end if
! Save coordinates for tallying purposes
p % last_xyz = p % coord(1) % xyz
! If particle has too many events, display warning and kill it
n_event = n_event + 1
if (n_event == MAX_EVENTS) then

38
src/trigger.F90
View file

@ -12,6 +12,7 @@ module trigger
use mesh_header, only: RegularMesh
use trigger_header, only: TriggerObject
use tally, only: TallyObject
use tally_filter, only: MeshFilter
implicit none
@ -165,7 +166,7 @@ contains
else
! Initialize bins, filter level
matching_bins(1:t % n_filters) = 0
matching_bins(1:size(t % filters)) = 0
FILTER_LOOP: do filter_index = 1, t % total_filter_bins
@ -265,7 +266,7 @@ contains
end if
end if
end do NUCLIDE_LOOP
if (t % n_filters == 0) exit FILTER_LOOP
if (size(t % filters) == 0) exit FILTER_LOOP
end do FILTER_LOOP
end if
end do TRIGGER_LOOP
@ -300,16 +301,19 @@ contains
! Get pointer to mesh
i_filter_mesh = t % find_filter(FILTER_MESH)
i_filter_surf = t % find_filter(FILTER_SURFACE)
m => meshes(t % filters(i_filter_mesh) % int_bins(1))
select type(filt => t % filters(i_filter_mesh) % obj)
type is (MeshFilter)
m => meshes(filt % mesh)
end select
! initialize bins array
matching_bins(1:t % n_filters) = 1
matching_bins(1:size(t % filters)) = 1
! determine how many energyin bins there are
i_filter_ein = t % find_filter(FILTER_ENERGYIN)
if (i_filter_ein > 0) then
print_ebin = .true.
n = t % filters(i_filter_ein) % n_bins
n = t % filters(i_filter_ein) % obj % n_bins
else
print_ebin = .false.
n = 1
@ -329,7 +333,7 @@ contains
mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_RIGHT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -341,7 +345,7 @@ contains
matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -356,7 +360,7 @@ contains
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_RIGHT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -368,7 +372,7 @@ contains
matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -383,7 +387,7 @@ contains
mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_FRONT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -396,7 +400,7 @@ contains
matching_bins(i_filter_surf) = OUT_FRONT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -411,7 +415,7 @@ contains
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_FRONT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -423,7 +427,7 @@ contains
matching_bins(i_filter_surf) = OUT_FRONT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -438,7 +442,7 @@ contains
mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.)
matching_bins(i_filter_surf) = IN_TOP
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -450,7 +454,7 @@ contains
matching_bins(i_filter_surf) = OUT_TOP
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -465,7 +469,7 @@ contains
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_TOP
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
@ -477,7 +481,7 @@ contains
matching_bins(i_filter_surf) = OUT_TOP
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
sum((matching_bins(1:size(t % filters)) - 1) * t % stride) + 1
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev

2
tests/test_score_current/results_true.dat
View file

@ -1 +1 @@
a9310752363eb059ff40f16ac9716b41ccab6ec6607d29f498069318745e485d18d784264304cc2586865bd58cef7587203cc22a1d485c58ddd63c14c0defdb9
bafab1921a12146abb2bb29603b52b9cc28a5a950a7a6bb1e3f012c05891c310fad643760d4f148b04d0fef3d1f3e141d146e3a278d81cc6fc8187c37717c5e7

2
tests/test_tallies/inputs_true.dat
View file

@ -1 +1 @@
ca47172a42f6c13b244a763c990cbe4811662708ee03307d810a4542ee34bb5db7cc29d66aea313dad95b9f38a4ff7943ded527cfd0c7c8825372fec40cfc0d0
930af242a043f2676a000dbc5a2db6b148edcb31ed8c87dbaa35a8efb37a3be8cff30cdf4dc03f9c5c7eb4021f7e4c3327e64681cdd8fd8722c95c69db850227

2
tests/test_tallies/results_true.dat
View file

@ -1 +1 @@
7d085e38f331083a8c3f7814eb6025f9457a1a8e3899e7d9a65af1ba92ea75c0182fcc605d4eadd3eb7edfc0949df688d39a3f45683e10cbdcfb6d7954f34129
a51db2a4efc681805f85968e04411dc33beee0532c202f5179b9a82880ab60a75e53fa9141c81045ea1d2842372f2d8da900326f09382ea61dd80a3c9b43bba1

7
tests/test_tallies/test_tallies.py
View file

@ -113,10 +113,11 @@ class TalliesTestHarness(PyAPITestHarness):
polar_tally4.estimator = 'tracklength'
universe_tally = Tally()
universe_tally.filters = [Filter(type='universe', bins=(1, 2, 3, 4))]
universe_tally.filters = [
Filter(type='universe', bins=(1, 2, 3, 4, 6, 8))]
universe_tally.scores = ['total']
cell_filter = Filter(type='cell', bins=(10, 21, 22, 23))
cell_filter = Filter(type='cell', bins=(10, 21, 22, 23, 60))
score_tallies = [Tally(), Tally(), Tally()]
for t in score_tallies:
t.filters = [cell_filter]
@ -128,7 +129,7 @@ class TalliesTestHarness(PyAPITestHarness):
score_tallies[1].estimator = 'analog'
score_tallies[2].estimator = 'collision'
cell_filter2 = Filter(type='cell', bins=(21, 22, 23, 27, 28, 29))
cell_filter2 = Filter(type='cell', bins=(21, 22, 23, 27, 28, 29, 60))
flux_tallies = [Tally() for i in range(4)]
for t in flux_tallies:
t.filters = [cell_filter2]