Get rid of a bunch of globals, add openmc_extend_meshes, revert to storing int

in mesh filter
This commit is contained in:
Paul Romano 2017-08-11 07:37:34 -05:00
parent c7985aa269
commit 01165aedad
14 changed files with 135 additions and 132 deletions

View file

@ -199,9 +199,6 @@ contains
n_source_points = 0
n_state_points = 0
n_tallies = 0
n_user_filters = 0
n_user_meshes = 0
n_user_tallies = 0
output_summary = .true.
output_tallies = .true.
particle_restart_run = .false.

View file

@ -50,14 +50,14 @@ contains
subroutine compute_xs()
use cmfd_header, only: cmfd_tallies
use constants, only: FILTER_MESH, FILTER_ENERGYIN, FILTER_ENERGYOUT, &
FILTER_SURFACE, OUT_LEFT, OUT_RIGHT, OUT_BACK, &
OUT_FRONT, OUT_BOTTOM, OUT_TOP, IN_LEFT, IN_RIGHT, &
IN_BACK, IN_FRONT, IN_BOTTOM, IN_TOP, CMFD_NOACCEL, &
ZERO, ONE, TINY_BIT
use error, only: fatal_error
use global, only: cmfd, n_cmfd_tallies, cmfd_tallies, meshes, &
filters, filter_matches
use global, only: cmfd, meshes, filters, filter_matches
use mesh_header, only: RegularMesh
use string, only: to_str
use tally_header, only: TallyObject
@ -75,7 +75,6 @@ contains
integer :: ital ! tally object index
integer :: ijk(3) ! indices for mesh cell
integer :: score_index ! index to pull from tally object
integer :: i_mesh ! index in meshes array
integer :: i_filt ! index in filters array
integer :: i_filter_mesh ! index for mesh filter
integer :: i_filter_ein ! index for incoming energy filter
@ -102,7 +101,7 @@ contains
i_filt = t % filter(t % find_filter(FILTER_MESH))
select type(filt => filters(i_filt) % obj)
type is (MeshFilter)
m => filt % mesh
m => meshes(filt % mesh)
end select
! Set mesh widths
@ -113,14 +112,14 @@ contains
cmfd % keff_bal = ZERO
! Begin loop around tallies
TAL: do ital = 1, n_cmfd_tallies
TAL: do ital = 1, size(cmfd_tallies)
! Associate tallies and mesh
t => cmfd_tallies(ital)
i_filt = t % filter(t % find_filter(FILTER_MESH))
select type(filt => filters(i_filt) % obj)
type is (MeshFilter)
m => filt % mesh
m => meshes(filt % mesh)
end select
! Check for energy filters

View file

@ -214,9 +214,10 @@ contains
subroutine cmfd_reweight(new_weights)
use algorithm, only: binary_search
use cmfd_header, only: cmfd_mesh
use constants, only: ZERO, ONE
use error, only: warning, fatal_error
use global, only: meshes, source_bank, work, n_user_meshes, cmfd
use global, only: source_bank, work, cmfd
use mesh_header, only: RegularMesh
use mesh, only: count_bank_sites
use message_passing
@ -234,14 +235,10 @@ contains
integer :: n_groups ! number of energy groups
logical :: outside ! any source sites outside mesh
logical :: in_mesh ! source site is inside mesh
type(RegularMesh), pointer :: m ! point to mesh
#ifdef MPI
integer :: mpi_err
#endif
! Associate pointer
m => meshes(n_user_meshes + 1)
! Get maximum of spatial and group indices
nx = cmfd % indices(1)
ny = cmfd % indices(2)
@ -265,8 +262,8 @@ contains
cmfd%weightfactors = ONE
! Count bank sites in mesh and reverse due to egrid structure
call count_bank_sites(m, source_bank, cmfd%sourcecounts, cmfd % egrid, &
sites_outside=outside, size_bank=work)
call count_bank_sites(cmfd_mesh, source_bank, cmfd%sourcecounts, &
cmfd % egrid, sites_outside=outside, size_bank=work)
cmfd % sourcecounts = cmfd%sourcecounts(ng:1:-1,:,:,:)
! Check for sites outside of the mesh
@ -295,7 +292,7 @@ contains
do i = 1, int(work,4)
! Determine spatial bin
call m % get_indices(source_bank(i) % xyz, ijk, in_mesh)
call cmfd_mesh % get_indices(source_bank(i) % xyz, ijk, in_mesh)
! Determine energy bin
n_groups = size(cmfd % egrid) - 1
@ -363,7 +360,7 @@ contains
subroutine cmfd_tally_reset()
use global, only: cmfd_tallies
use cmfd_header, only: cmfd_tallies
use output, only: write_message
integer :: i ! loop counter

View file

@ -1,6 +1,8 @@
module cmfd_header
use constants, only: CMFD_NOACCEL, ZERO, ONE
use mesh_header, only: RegularMesh
use tally_header, only: TallyObject
implicit none
private
@ -88,6 +90,11 @@ module cmfd_header
end type cmfd_type
type(RegularMesh), public, pointer :: cmfd_mesh => null()
! Pointers for different tallies
type(TallyObject), public, pointer :: cmfd_tallies(:) => null()
contains
!==============================================================================

View file

@ -248,9 +248,10 @@ contains
subroutine create_cmfd_tally(root)
use cmfd_header, only: cmfd_mesh, cmfd_tallies
use constants, only: MAX_LINE_LEN
use error, only: fatal_error, warning
use mesh_header, only: RegularMesh
use mesh_header, only: RegularMesh, openmc_extend_meshes
use string
use tally_header, only: TallyObject, openmc_extend_tallies
use tally_filter_header
@ -265,6 +266,7 @@ contains
integer :: ng ! number of energy groups (default 1)
integer :: n_filter ! number of filters
integer :: i_start, i_end
integer :: i_filt_start, i_filt_end
integer(C_INT) :: err
integer :: i_filt ! index in filters array
integer :: iarray3(3) ! temp integer array
@ -273,16 +275,14 @@ contains
type(RegularMesh), pointer :: m
type(XMLNode) :: node_mesh
! Set global variables if they are 0 (this can happen if there is no tally
! file)
if (n_meshes == 0) n_meshes = n_user_meshes + n_cmfd_meshes
err = openmc_extend_meshes(1, i_start)
! Allocate mesh
if (.not. allocated(meshes)) allocate(meshes(n_meshes))
m => meshes(n_user_meshes+1)
cmfd_mesh => meshes(i_start)
m => meshes(i_start)
! Set mesh id
m % id = n_user_meshes + 1
m % id = i_start
! Set mesh type to rectangular
m % type = LATTICE_RECT
@ -376,23 +376,23 @@ contains
m % volume_frac = ONE/real(product(m % dimension),8)
! Add mesh to dictionary
call mesh_dict % add_key(m % id, n_user_meshes + 1)
call mesh_dict % add_key(m % id, i_start)
! Determine number of filters
energy_filters = check_for_node(node_mesh, "energy")
n_cmfd_filters = merge(5, 3, energy_filters)
n = merge(5, 3, energy_filters)
! Extend filters array so we can add CMFD filters
err = openmc_extend_filters(n_cmfd_filters, i_start, i_end)
err = openmc_extend_filters(n, i_filt_start, i_filt_end)
! Set up mesh filter
i_filt = i_start
i_filt = i_filt_start
allocate(MeshFilter :: filters(i_filt) % obj)
select type (filt => filters(i_filt) % obj)
type is (MeshFilter)
filt % id = i_filt
filt % n_bins = product(m % dimension)
filt % mesh => meshes(n_user_meshes + 1)
filt % mesh = i_start
! Add filter to dictionary
call filter_dict % add_key(filt % id, i_filt)
end select
@ -437,7 +437,7 @@ contains
! 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)
filt % mesh => meshes(n_user_meshes + 1)
filt % mesh = i_start
! Add filter to dictionary
call filter_dict % add_key(filt % id, i_filt)
end select
@ -464,11 +464,11 @@ contains
end select
! Allocate tallies
err = openmc_extend_tallies(n_cmfd_tallies, i_start, i_end)
err = openmc_extend_tallies(3, i_start, i_end)
cmfd_tallies => tallies(i_start:i_end)
! Begin loop around tallies
do i = 1, n_cmfd_tallies
do i = 1, size(cmfd_tallies)
! Point t to tally variable
t => cmfd_tallies(i)
@ -510,9 +510,9 @@ contains
! Allocate and set filters
allocate(t % filter(n_filter))
t % filter(1) = n_user_filters + 1
t % filter(1) = i_filt_start
if (energy_filters) then
t % filter(2) = n_user_filters + 2
t % filter(2) = i_filt_start + 1
end if
! Allocate scoring bins
@ -550,10 +550,10 @@ contains
! Allocate and set indices in filters array
allocate(t % filter(n_filter))
t % filter(1) = n_user_filters + 1
t % filter(1) = i_filt_start
if (energy_filters) then
t % filter(2) = n_user_filters + 2
t % filter(3) = n_user_filters + 3
t % filter(2) = i_filt_start + 1
t % filter(3) = i_filt_start + 2
end if
! Allocate macro reactions
@ -583,10 +583,10 @@ contains
! Allocate and set filters
allocate(t % filter(n_filter))
t % filter(1) = n_user_filters + n_cmfd_filters - 1
t % filter(n_filter) = n_user_filters + n_cmfd_filters
t % filter(1) = i_filt_end - 1
t % filter(n_filter) = i_filt_end
if (energy_filters) then
t % filter(2) = n_user_filters + 2
t % filter(2) = i_filt_start + 1
end if
! Allocate macro reactions

View file

@ -71,10 +71,6 @@ module global
! ============================================================================
! TALLY-RELATED VARIABLES
! Pointers for different tallies
type(TallyObject), pointer :: user_tallies(:) => null()
type(TallyObject), pointer :: cmfd_tallies(:) => null()
! Active tally lists
type(VectorInt) :: active_analog_tallies
type(VectorInt) :: active_tracklength_tallies
@ -83,10 +79,6 @@ module global
type(VectorInt) :: active_tallies
type(VectorInt) :: active_surface_tallies
integer :: n_user_meshes = 0 ! # of structured user meshes
integer :: n_user_filters = 0 ! # of user filters
integer :: n_user_tallies = 0 ! # of user tallies
! Normalization for statistics
integer :: n_realizations = 0 ! # of independent realizations
real(8) :: total_weight ! total starting particle weight in realization
@ -264,11 +256,6 @@ module global
! Flag to activate neutronic feedback via source weights
logical :: cmfd_feedback = .false.
! User-defined tally information
integer :: n_cmfd_meshes = 1 ! # of structured meshes
integer :: n_cmfd_filters = 0 ! # of filters
integer :: n_cmfd_tallies = 3 ! # of user-defined tallies
! Adjoint method type
character(len=10) :: cmfd_adjoint_type = 'physical'

View file

@ -4,6 +4,7 @@ module input_xml
use algorithm, only: find
use cmfd_input, only: configure_cmfd
use cmfd_header, only: cmfd_mesh
use constants
use dict_header, only: DictIntInt, DictCharInt, ElemKeyValueCI
use distribution_multivariate
@ -2652,6 +2653,7 @@ contains
integer :: trig_ind ! index of triggers array for each tally
integer :: user_trig_ind ! index of user-specified triggers for each tally
integer :: i_start, i_end
integer :: i_filt_start, i_filt_end
integer(C_INT) :: err
real(8) :: threshold ! trigger convergence threshold
integer :: n_order ! moment order requested
@ -2723,37 +2725,6 @@ contains
! Get pointer list to XML <tally>
call get_node_list(root, "tally", node_tal_list)
! Check for user meshes
n_user_meshes = size(node_mesh_list)
if (cmfd_run) then
n_meshes = n_user_meshes + n_cmfd_meshes
else
n_meshes = n_user_meshes
end if
! Allocate mesh array
if (n_meshes > 0) allocate(meshes(n_meshes))
! Check for user filters
n_user_filters = size(node_filt_list)
! Allocate filters array
if (n_user_filters > 0) then
err = openmc_extend_filters(n_user_filters)
end if
! Check for user tallies
n_user_tallies = size(node_tal_list)
if (n_user_tallies == 0) then
if (master) call warning("No tallies present in tallies.xml file!")
end if
! Allocate user tallies
if (n_user_tallies > 0 .and. run_mode /= MODE_PLOTTING) then
err = openmc_extend_tallies(n_user_tallies, i_start, i_end)
user_tallies => tallies(i_start:i_end)
end if
! Check for <assume_separate> setting
if (check_for_node(root, "assume_separate")) then
call get_node_value(root, "assume_separate", assume_separate)
@ -2762,8 +2733,14 @@ contains
! ==========================================================================
! READ MESH DATA
do i = 1, n_user_meshes
m => meshes(i)
! Check for user meshes and allocate
n = size(node_mesh_list)
if (n > 0) then
err = openmc_extend_meshes(n, i_start, i_end)
end if
do i = 1, n
m => meshes(i_start + i - 1)
! Get pointer to mesh node
node_mesh = node_mesh_list(i)
@ -2979,8 +2956,14 @@ contains
! ==========================================================================
! READ FILTER DATA
READ_FILTERS: do i = 1, n_user_filters
f => filters(i)
! Check for user filters and allocate
n = size(node_filt_list)
if (n > 0) then
err = openmc_extend_filters(n, i_start, i_end)
end if
READ_FILTERS: do i = 1, n
f => filters(i_start + i - 1)
! Get pointer to filter xml node
node_filt = node_filt_list(i)
@ -3123,7 +3106,7 @@ contains
filt % n_bins = product(m % dimension)
! Store the index of the mesh
filt % mesh => meshes(i_mesh)
filt % mesh = i_mesh
end select
case ('energy')
@ -3342,9 +3325,20 @@ contains
! ==========================================================================
! READ TALLY DATA
READ_TALLIES: do i = 1, n_user_tallies
! Check for user tallies
n = size(node_tal_list)
if (n == 0) then
if (master) call warning("No tallies present in tallies.xml file!")
end if
! Allocate user tallies
if (n > 0 .and. run_mode /= MODE_PLOTTING) then
err = openmc_extend_tallies(n, i_start, i_end)
end if
READ_TALLIES: do i = 1, n
! Get pointer to tally
t => tallies(i)
t => tallies(i_start + i - 1)
! Get pointer to tally xml node
node_tal = node_tal_list(i)
@ -3865,7 +3859,8 @@ contains
! Get pointer to mesh
select type(filt => filters(i_filter_mesh) % obj)
type is (MeshFilter)
m => filt % mesh
i_mesh = filt % mesh
m => meshes(filt % mesh)
end select
! Copy filter indices to temporary array
@ -3879,13 +3874,10 @@ contains
! Extend the filters array so we can add a surface filter and
! mesh filter
err = openmc_extend_filters(2, i_start, i_end)
! Increment number of user filters
n_user_filters = n_user_filters + 2
err = openmc_extend_filters(2, i_filt_start, i_filt_end)
! Get index of the new mesh filter
i_filt = i_start
i_filt = i_filt_start
! Duplicate the mesh filter since other tallies might use this
! filter and we need to change the dimension
@ -3893,7 +3885,7 @@ contains
select type(filt => filters(i_filt) % obj)
type is (MeshFilter)
filt % id = i_filt
filt % mesh => m
filt % mesh = i_mesh
! We need to increase the dimension by one since we also need
! currents coming into and out of the boundary mesh cells.
@ -3905,7 +3897,7 @@ contains
t % filter(t % find_filter(FILTER_MESH)) = i_filt
! Get index of the new surface filter
i_filt = i_end
i_filt = i_filt_end
! Add surface filter
allocate(SurfaceFilter :: filters(i_filt) % obj)
@ -4683,11 +4675,7 @@ contains
&meshlines on plot " // trim(to_str(pl % id)))
end if
select type(filt => filters(cmfd_tallies(1) % &
filter(cmfd_tallies(1) % find_filter(FILTER_MESH))) % obj)
type is (MeshFilter)
pl % meshlines_mesh => filt % mesh
end select
pl % meshlines_mesh => cmfd_mesh
case ('entropy')

View file

@ -11,7 +11,7 @@ module mesh_header
implicit none
private
public :: RegularMesh, meshes, n_meshes, mesh_dict
public :: openmc_extend_meshes
!===============================================================================
! STRUCTUREDMESH represents a tessellation of n-dimensional Euclidean space by
@ -36,12 +36,12 @@ module mesh_header
procedure :: to_hdf5 => regular_to_hdf5
end type RegularMesh
integer(C_INT32_T), bind(C) :: n_meshes = 0 ! # of structured meshes
integer(C_INT32_T), public, bind(C) :: n_meshes = 0 ! # of structured meshes
type(RegularMesh), allocatable, target :: meshes(:)
type(RegularMesh), public, allocatable, target :: meshes(:)
! Dictionary that maps user IDs to indices in 'meshes'
type(DictIntInt) :: mesh_dict
type(DictIntInt), public :: mesh_dict
contains
@ -164,8 +164,8 @@ contains
pure function regular_intersects(this, xyz0, xyz1) result(intersects)
class(RegularMesh), intent(in) :: this
real(8), intent(in) :: xyz0(1)
real(8), intent(in) :: xyz1(1)
real(8), intent(in) :: xyz0(:)
real(8), intent(in) :: xyz1(:)
logical :: intersects
select case(this % n_dimension)
@ -415,4 +415,39 @@ contains
call close_group(mesh_group)
end subroutine regular_to_hdf5
!===============================================================================
! C API FUNCTIONS
!===============================================================================
function openmc_extend_meshes(n, index_start, index_end) result(err) bind(C)
! Extend the meshes array by n elements
integer(C_INT32_T), value, intent(in) :: n
integer(C_INT32_T), optional, intent(out) :: index_start
integer(C_INT32_T), optional, intent(out) :: index_end
integer(C_INT) :: err
type(RegularMesh), allocatable :: temp(:) ! temporary meshes array
if (n_meshes == 0) then
! Allocate meshes array
allocate(meshes(n))
else
! Allocate meshes array with increased size
allocate(temp(n_meshes + n))
! Copy original meshes to temporary array
temp(1:n_meshes) = meshes
! Move allocation from temporary array
call move_alloc(FROM=temp, TO=meshes)
end if
! Return indices in meshes array
if (present(index_start)) index_start = n_meshes + 1
if (present(index_end)) index_end = n_meshes + n
n_meshes = n_meshes + n
err = 0
end function openmc_extend_meshes
end module mesh_header

View file

@ -1023,7 +1023,7 @@ contains
i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
select type(filt => filters(i_filter_mesh) % obj)
type is (MeshFilter)
m => filt % mesh
m => meshes(filt % mesh)
end select
! Get surface filter index and stride

View file

@ -2223,7 +2223,6 @@ contains
integer :: k ! loop index for nuclide bins
integer :: filter_index ! single index for single bin
integer :: i_nuclide ! index in nuclides array
integer :: matching_bin ! next valid filter bin
real(8) :: filter_weight ! combined weight of all filters
logical :: finished ! found all valid bin combinations
type(TallyObject), pointer :: t
@ -2368,7 +2367,6 @@ contains
integer :: k ! loop index for nuclide bins
integer :: filter_index ! single index for single bin
integer :: i_nuclide ! index in nuclides array
integer :: matching_bin ! next valid filter bin
real(8) :: filter_weight ! combined weight of all filters
real(8) :: atom_density
logical :: finished ! found all valid bin combinations
@ -2744,7 +2742,6 @@ contains
integer :: k ! loop index for nuclide bins
integer :: filter_index ! single index for single bin
integer :: i_nuclide ! index in nuclides array (from bins)
integer :: matching_bin ! next valid filter bin
real(8) :: flux ! tracklength estimate of flux
real(8) :: atom_density ! atom density of single nuclide in atom/b-cm
real(8) :: filter_weight ! combined weight of all filters
@ -2905,7 +2902,6 @@ contains
integer :: k ! loop index for nuclide bins
integer :: filter_index ! single index for single bin
integer :: i_nuclide ! index in nuclides array (from bins)
integer :: matching_bin ! next valid filter bin
real(8) :: flux ! collision estimate of flux
real(8) :: atom_density ! atom density of single nuclide
! in atom/b-cm
@ -3192,7 +3188,6 @@ contains
integer :: i
integer :: i_tally
integer :: i_filt
integer :: j, k ! loop indices
integer :: n_dim ! num dimensions of the mesh
integer :: d1 ! dimension index
@ -3211,7 +3206,6 @@ contains
real(8) :: xyz_cross(3) ! coordinates of bounding surfaces
real(8) :: d(3) ! distance to each bounding surface
real(8) :: distance ! actual distance traveled
real(8) :: filt_score ! score applied by filters
integer :: matching_bin ! next valid filter bin
logical :: start_in_mesh ! particle's starting xyz in mesh?
logical :: end_in_mesh ! particle's ending xyz in mesh?
@ -3249,7 +3243,7 @@ contains
! Get pointer to mesh
select type(filt => filters(i_filter_mesh) % obj)
type is (MeshFilter)
m => filt % mesh
m => meshes(filt % mesh)
end select
n_dim = m % n_dimension

View file

@ -188,7 +188,7 @@ contains
integer, intent(in) :: estimator
type(TallyFilterMatch), intent(inout) :: match
integer :: i, start
integer :: i
! Iterate over coordinate levels to see which universes match
do i = 1, p % n_coord

View file

@ -4,7 +4,7 @@ module tally_filter_mesh
use constants
use error, only: warning
use mesh_header, only: RegularMesh
use mesh_header, only: RegularMesh, meshes
use hdf5_interface
use particle_header, only: Particle
use string, only: to_str
@ -20,7 +20,7 @@ module tally_filter_mesh
!===============================================================================
type, public, extends(TallyFilter) :: MeshFilter
type(RegularMesh), pointer :: mesh => null()
integer :: mesh
contains
procedure :: get_all_bins => get_all_bins_mesh
procedure :: to_statepoint => to_statepoint_mesh
@ -60,7 +60,7 @@ contains
weight = ERROR_REAL
! Get a pointer to the mesh.
m => this % mesh
m => meshes(this % mesh)
n = m % n_dimension
if (estimator /= ESTIMATOR_TRACKLENGTH) then
@ -215,7 +215,7 @@ contains
call write_dataset(filter_group, "type", "mesh")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", this % mesh % id)
call write_dataset(filter_group, "bins", meshes(this % mesh) % id)
end subroutine to_statepoint_mesh
function text_label_mesh(this, bin) result(label)
@ -225,7 +225,7 @@ contains
integer, allocatable :: ijk(:)
associate (m => this % mesh)
associate (m => meshes(this % mesh))
allocate(ijk(m % n_dimension))
call m % get_indices_from_bin(bin, ijk)
if (m % n_dimension == 1) then

View file

@ -213,7 +213,6 @@ contains
subroutine tally_setup_arrays(this)
class(TallyObject), intent(inout) :: this
integer :: i ! loop index for tallies
integer :: j ! loop index for filters
integer :: n ! temporary stride
integer :: i_filt ! filter index
@ -274,8 +273,8 @@ contains
function openmc_extend_tallies(n, index_start, index_end) result(err) bind(C)
! Extend the tallies array by n elements
integer(C_INT32_T), value, intent(in) :: n
integer(C_INT32_T), intent(out) :: index_start
integer(C_INT32_T), intent(out) :: index_end
integer(C_INT32_T), optional, intent(out) :: index_start
integer(C_INT32_T), optional, intent(out) :: index_end
integer(C_INT) :: err
type(TallyObject), allocatable :: temp(:) ! temporary tallies array
@ -295,9 +294,9 @@ contains
end if
! Return indices in tallies array
index_start = n_tallies + 1
index_end = n_tallies + n
n_tallies = index_end
if (present(index_start)) index_start = n_tallies + 1
if (present(index_end)) index_end = n_tallies + n
n_tallies = n_tallies + n
err = 0
end function openmc_extend_tallies

View file

@ -315,7 +315,7 @@ contains
i_filter_surf = t % filter(t % find_filter(FILTER_SURFACE))
select type(filt => filters(i_filter_mesh) % obj)
type is (MeshFilter)
m => filt % mesh
m => meshes(filt % mesh)
end select
! initialize bins array