Store entropy and UFS meshes in global meshes array

This commit is contained in:
Paul Romano 2017-08-21 11:16:56 -05:00
parent ee70dfdeca
commit 569a532af6
6 changed files with 107 additions and 141 deletions

View file

@ -300,39 +300,9 @@ contains
integer :: ent_idx ! entropy index
integer :: i, j, k ! index for bank sites
integer :: n ! # of boxes in each dimension
logical :: sites_outside ! were there sites outside entropy box?
type(RegularMesh), pointer :: m
! Get pointer to entropy mesh
m => entropy_mesh
! On the first pass through this subroutine, we need to determine how big
! the entropy mesh should be in each direction and then allocate a
! three-dimensional array to store the fraction of source sites in each mesh
! box
if (.not. allocated(entropy_p)) then
if (.not. allocated(m % dimension)) then
! If the user did not specify how many mesh cells are to be used in
! each direction, we automatically determine an appropriate number of
! cells
n = ceiling((n_particles/20)**(ONE/THREE))
! copy dimensions
m % n_dimension = 3
allocate(m % dimension(3))
m % dimension = n
! determine width
m % width = (m % upper_right - m % lower_left) / m % dimension
end if
! allocate p
allocate(entropy_p(1, m % dimension(1), m % dimension(2), &
m % dimension(3)))
end if
associate (m => meshes(index_entropy_mesh))
! count number of fission sites over mesh
call count_bank_sites(m, fission_bank, entropy_p, &
@ -362,6 +332,7 @@ contains
end do
end if
end associate
end subroutine shannon_entropy
!===============================================================================
@ -621,16 +592,18 @@ contains
integer :: mpi_err ! MPI error code
#endif
associate (m => meshes(index_ufs_mesh))
if (current_batch == 1 .and. current_gen == 1) then
! On the first generation, just assume that the source is already evenly
! distributed so that effectively the production of fission sites is not
! biased
source_frac = ufs_mesh % volume_frac
source_frac = m % volume_frac
else
! count number of source sites in each ufs mesh cell
call count_bank_sites(ufs_mesh, source_bank, source_frac, &
call count_bank_sites(m, source_bank, source_frac, &
sites_outside=sites_outside, size_bank=work)
! Check for sites outside of the mesh
@ -640,7 +613,7 @@ contains
#ifdef MPI
! Send source fraction to all processors
n = product(ufs_mesh % dimension)
n = product(m % dimension)
call MPI_BCAST(source_frac, n, MPI_REAL8, 0, mpi_intracomm, mpi_err)
#endif
@ -654,6 +627,8 @@ contains
source_bank % wgt = source_bank % wgt * n_particles / total
end if
end associate
end subroutine count_source_for_ufs
#ifdef _OPENMP

View file

@ -202,25 +202,8 @@ contains
call statepoint_batch % clear()
call sourcepoint_batch % clear()
! Deallocate entropy mesh
if (associated(entropy_mesh)) then
if (allocated(entropy_mesh % lower_left)) &
deallocate(entropy_mesh % lower_left)
if (allocated(entropy_mesh % upper_right)) &
deallocate(entropy_mesh % upper_right)
if (allocated(entropy_mesh % width)) deallocate(entropy_mesh % width)
deallocate(entropy_mesh)
end if
! Deallocate ufs
if (allocated(source_frac)) deallocate(source_frac)
if (associated(ufs_mesh)) then
if (allocated(ufs_mesh % lower_left)) deallocate(ufs_mesh % lower_left)
if (allocated(ufs_mesh % upper_right)) &
deallocate(ufs_mesh % upper_right)
if (allocated(ufs_mesh % width)) deallocate(ufs_mesh % width)
deallocate(ufs_mesh)
end if
end subroutine free_memory

View file

@ -144,6 +144,7 @@ contains
integer :: n
integer :: temp_int
integer :: temp_int_array3(3)
integer(C_INT) :: err
integer, allocatable :: temp_int_array(:)
real(8), allocatable :: temp_real(:)
integer :: n_tracks
@ -700,47 +701,54 @@ contains
&corner of Shannon entropy mesh.")
end if
! Allocate mesh object and coordinates on mesh
allocate(entropy_mesh)
allocate(entropy_mesh % lower_left(3))
allocate(entropy_mesh % upper_right(3))
allocate(entropy_mesh % width(3))
err = openmc_extend_meshes(1, index_entropy_mesh)
! Copy values
call get_node_array(node_entropy, "lower_left", &
entropy_mesh % lower_left)
call get_node_array(node_entropy, "upper_right", &
entropy_mesh % upper_right)
associate (m => meshes(index_entropy_mesh))
! Allocate mesh object and coordinates on mesh
allocate(m % lower_left(3))
allocate(m % upper_right(3))
allocate(m % width(3))
! Check on values provided
if (.not. all(entropy_mesh % upper_right > entropy_mesh % lower_left)) &
&then
call fatal_error("Upper-right coordinate must be greater than &
&lower-left coordinate for Shannon entropy mesh.")
end if
! Copy values
call get_node_array(node_entropy, "lower_left", m % lower_left)
call get_node_array(node_entropy, "upper_right", m % upper_right)
! Check if dimensions were specified -- if not, they will be calculated
! automatically upon first entry into shannon_entropy
if (check_for_node(node_entropy, "dimension")) then
! If so, make sure proper number of values were given
if (node_word_count(node_entropy, "dimension") /= 3) then
call fatal_error("Dimension of entropy mesh must be given as three &
&integers.")
! Check on values provided
if (.not. all(m % upper_right > m % lower_left)) &
&then
call fatal_error("Upper-right coordinate must be greater than &
&lower-left coordinate for Shannon entropy mesh.")
end if
! Allocate dimensions
entropy_mesh % n_dimension = 3
allocate(entropy_mesh % dimension(3))
m % n_dimension = 3
allocate(m % dimension(3))
! Copy dimensions
call get_node_array(node_entropy, "dimension", entropy_mesh % dimension)
! Check if dimensions were specified -- if not, they will be calculated
! automatically upon first entry into shannon_entropy
if (check_for_node(node_entropy, "dimension")) then
! If so, make sure proper number of values were given
if (node_word_count(node_entropy, "dimension") /= 3) then
call fatal_error("Dimension of entropy mesh must be given as three &
&integers.")
end if
! Copy dimensions
call get_node_array(node_entropy, "dimension", m % dimension)
else
! If the user did not specify how many mesh cells are to be used in
! each direction, we automatically determine an appropriate number of
! cells
m % dimension = ceiling((n_particles/20)**(ONE/THREE))
end if
! Calculate width
entropy_mesh % width = (entropy_mesh % upper_right - &
entropy_mesh % lower_left) / entropy_mesh % dimension
m % width = (m % upper_right - m % lower_left) / m % dimension
end if
! Allocate space for storing number of fission sites in each mesh cell
allocate(entropy_p(1, m % dimension(1), m % dimension(2), &
m % dimension(3)))
end associate
! Turn on Shannon entropy calculation
entropy_on = .true.
@ -764,42 +772,44 @@ contains
&integers.")
end if
err = openmc_extend_meshes(1, index_ufs_mesh)
! Allocate mesh object and coordinates on mesh
allocate(ufs_mesh)
allocate(ufs_mesh % lower_left(3))
allocate(ufs_mesh % upper_right(3))
allocate(ufs_mesh % width(3))
associate (m => meshes(index_ufs_mesh))
allocate(m % lower_left(3))
allocate(m % upper_right(3))
allocate(m % width(3))
! Allocate dimensions
ufs_mesh % n_dimension = 3
allocate(ufs_mesh % dimension(3))
! Allocate dimensions
m % n_dimension = 3
allocate(m % dimension(3))
! Copy dimensions
call get_node_array(node_ufs, "dimension", ufs_mesh % dimension)
! Copy dimensions
call get_node_array(node_ufs, "dimension", m % dimension)
! Copy values
call get_node_array(node_ufs, "lower_left", ufs_mesh % lower_left)
call get_node_array(node_ufs, "upper_right", ufs_mesh % upper_right)
! Copy values
call get_node_array(node_ufs, "lower_left", m % lower_left)
call get_node_array(node_ufs, "upper_right", m % upper_right)
! Check on values provided
if (.not. all(ufs_mesh % upper_right > ufs_mesh % lower_left)) then
call fatal_error("Upper-right coordinate must be greater than &
&lower-left coordinate for UFS mesh.")
end if
! Check on values provided
if (.not. all(m % upper_right > m % lower_left)) then
call fatal_error("Upper-right coordinate must be greater than &
&lower-left coordinate for UFS mesh.")
end if
! Calculate width
ufs_mesh % width = (ufs_mesh % upper_right - &
ufs_mesh % lower_left) / ufs_mesh % dimension
! Calculate width
m % width = (m % upper_right - m % lower_left) / m % dimension
! Calculate volume fraction of each cell
ufs_mesh % volume_frac = ONE/real(product(ufs_mesh % dimension),8)
! Calculate volume fraction of each cell
m % volume_frac = ONE/real(product(m % dimension),8)
! Allocate source_frac
allocate(source_frac(1, m % dimension(1), m % dimension(2), &
m % dimension(3)))
end associate
! Turn on uniform fission source weighting
ufs = .true.
! Allocate source_frac
allocate(source_frac(1, ufs_mesh % dimension(1), &
ufs_mesh % dimension(2), ufs_mesh % dimension(3)))
end if
! Check if the user has specified to write state points
@ -4566,12 +4576,12 @@ contains
select case (trim(meshtype))
case ('ufs')
if (.not. associated(ufs_mesh)) then
if (index_ufs_mesh < 0) then
call fatal_error("No UFS mesh for meshlines on plot " &
// trim(to_str(pl % id)))
end if
pl % meshlines_mesh => ufs_mesh
pl % meshlines_mesh => meshes(index_ufs_mesh)
case ('cmfd')
@ -4584,17 +4594,12 @@ contains
case ('entropy')
if (.not. associated(entropy_mesh)) then
if (index_entropy_mesh < 0) then
call fatal_error("No entropy mesh for meshlines on plot " &
// trim(to_str(pl % id)))
end if
if (.not. allocated(entropy_mesh % dimension)) then
call fatal_error("No dimension specified on entropy mesh &
&for meshlines on plot " // trim(to_str(pl % id)))
end if
pl % meshlines_mesh => entropy_mesh
pl % meshlines_mesh => meshes(index_entropy_mesh)
case ('tally')

View file

@ -20,7 +20,7 @@ contains
subroutine count_bank_sites(m, bank_array, cnt, energies, size_bank, &
sites_outside)
type(RegularMesh), pointer :: m ! mesh to count sites
type(RegularMesh), intent(in) :: m ! mesh to count sites
type(Bank), intent(in) :: bank_array(:) ! fission or source bank
real(8), intent(out) :: cnt(:,:,:,:) ! weight of sites in each
! cell and energy group

View file

@ -1088,18 +1088,20 @@ contains
! the expected number of fission sites produced
if (ufs) then
! Determine indices on ufs mesh for current location
call ufs_mesh % get_indices(p % coord(1) % xyz, ijk, in_mesh)
if (.not. in_mesh) then
call write_particle_restart(p)
call fatal_error("Source site outside UFS mesh!")
end if
associate (m => meshes(index_ufs_mesh))
! Determine indices on ufs mesh for current location
call m % get_indices(p % coord(1) % xyz, ijk, in_mesh)
if (.not. in_mesh) then
call write_particle_restart(p)
call fatal_error("Source site outside UFS mesh!")
end if
if (source_frac(1,ijk(1),ijk(2),ijk(3)) /= ZERO) then
weight = ufs_mesh % volume_frac / source_frac(1,ijk(1),ijk(2),ijk(3))
else
weight = ONE
end if
if (source_frac(1,ijk(1),ijk(2),ijk(3)) /= ZERO) then
weight = m % volume_frac / source_frac(1,ijk(1),ijk(2),ijk(3))
else
weight = ONE
end if
end associate
else
weight = ONE
end if

View file

@ -193,20 +193,21 @@ contains
! the expected number of fission sites produced
if (ufs) then
associate (m => meshes(index_ufs_mesh))
! Determine indices on ufs mesh for current location
call m % get_indices(p % coord(1) % xyz, ijk, in_mesh)
! Determine indices on ufs mesh for current location
call ufs_mesh % get_indices(p % coord(1) % xyz, ijk, in_mesh)
if (.not. in_mesh) then
call write_particle_restart(p)
call fatal_error("Source site outside UFS mesh!")
end if
if (.not. in_mesh) then
call write_particle_restart(p)
call fatal_error("Source site outside UFS mesh!")
end if
if (source_frac(1,ijk(1),ijk(2),ijk(3)) /= ZERO) then
weight = ufs_mesh % volume_frac / source_frac(1,ijk(1),ijk(2),ijk(3))
else
weight = ONE
end if
if (source_frac(1,ijk(1),ijk(2),ijk(3)) /= ZERO) then
weight = m % volume_frac / source_frac(1,ijk(1),ijk(2),ijk(3))
else
weight = ONE
end if
end associate
else
weight = ONE
end if