diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index f2cb4f5e12..587fb06cca 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -148,7 +148,10 @@ main.o: string.o main.o: tally.o main.o: timing.o +mesh.o: constants.o +mesh.o: global.o mesh.o: mesh_header.o +mesh.o: particle_header.o output.o: constants.o output.o: datatypes.o diff --git a/src/global.F90 b/src/global.F90 index 5fadb2e1af..752e0db3a7 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -152,8 +152,8 @@ module global ! Shannon entropy logical :: entropy_on = .false. - real(8) :: entropy ! value of shannon entropy - real(8), allocatable :: entropy_p(:) ! fraction of source sites in each cell + real(8) :: entropy ! value of shannon entropy + real(8), allocatable :: entropy_p(:,:,:) ! % of source sites in each cell type(StructuredMesh), pointer :: entropy_mesh ! ============================================================================ diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 6d3ae131e8..97a7b9f86d 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -4,7 +4,7 @@ module intercycle use error, only: fatal_error, warning use global - use mesh, only: get_mesh_bin + use mesh, only: count_fission_sites use mesh_header, only: StructuredMesh use output, only: write_message use random_lcg, only: prn, set_particle_seed, prn_skip @@ -283,12 +283,10 @@ contains subroutine shannon_entropy() - integer :: i ! index for bank sites + integer :: i, j, k ! index for bank sites integer :: n ! # of boxes in each dimension - integer :: bin ! index in entropy_p - integer(8) :: total_bank ! total # of fission bank sites - integer, save :: n_box ! total # of boxes on mesh - logical :: outside_box ! were there sites outside entropy box? + real(8) :: total ! total weight of fission bank sites + logical :: sites_outside ! were there sites outside entropy box? type(StructuredMesh), pointer :: m => null() ! Get pointer to entropy mesh @@ -312,69 +310,39 @@ contains m % dimension = n end if - ! Determine total number of mesh boxes - n_box = product(m % dimension) - ! allocate and determine width allocate(m % width(3)) m % width = (m % upper_right - m % lower_left) / m % dimension ! allocate p - allocate(entropy_p(n_box)) + allocate(entropy_p(m % dimension(1), m % dimension(2), & + m % dimension(3))) end if - ! initialize p - entropy_p = ZERO - outside_box = .false. - - ! loop over fission sites and count how many are in each mesh box - FISSION_SITES: do i = 1, int(n_bank,4) - ! determine scoring bin for entropy mesh - call get_mesh_bin(m, fission_bank(i) % xyz, bin) - - ! if outside mesh, skip particle - if (bin == NO_BIN_FOUND) then - outside_box = .true. - cycle - end if - - ! add to appropriate mesh box - entropy_p(bin) = entropy_p(bin) + 1 - end do FISSION_SITES + ! count number of fission sites over mesh + call count_fission_sites(m, entropy_p, total, sites_outside) ! display warning message if there were sites outside entropy box - if (outside_box) then + if (sites_outside) then message = "Fission source site(s) outside of entropy box." call warning() end if -#ifdef MPI - ! collect values from all processors - if (master) then - call MPI_REDUCE(MPI_IN_PLACE, entropy_p, n_box, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) - else - call MPI_REDUCE(entropy_p, entropy_p, n_box, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) - end if - - ! determine total number of bank sites - call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, & - MPI_COMM_WORLD, mpi_err) -#else - total_bank = n_bank -#endif - ! sum values to obtain shannon entropy if (master) then - ! Normalize to number of bank sites - entropy_p = entropy_p / total_bank + ! Normalize to total weight of bank sites + entropy_p = entropy_p / total entropy = 0 - do i = 1, n_box - if (entropy_p(i) > 0) then - entropy = entropy - entropy_p(i) * log(entropy_p(i))/log(2.0) - end if + do i = 1, m % dimension(1) + do j = 1, m % dimension(2) + do k = 1, m % dimension(3) + if (entropy_p(i,j,k) > ZERO) then + entropy = entropy - entropy_p(i,j,k) * & + log(entropy_p(i,j,k))/log(2.0) + end if + end do + end do end do end if diff --git a/src/mesh.F90 b/src/mesh.F90 index 22c2f28c86..756525327c 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -1,9 +1,14 @@ module mesh use constants + use global use mesh_header use particle_header, only: Particle +#ifdef MPI + use mpi +#endif + implicit none contains @@ -137,4 +142,79 @@ contains end subroutine bin_to_mesh_indices +!=============================================================================== +! COUNT_FISSION_SITES determines the number of fission bank sites in each cell +! of a given mesh. This can be used for a variety of purposes (Shannon entropy, +! CMFD, uniform fission source weighting) +!=============================================================================== + + subroutine count_fission_sites(m, count, total, sites_outside) + + type(StructuredMesh), pointer :: m ! mesh to count sites + real(8), intent(out) :: count(:,:,:) ! weight of sites in each cell + real(8), intent(out) :: total ! total weight of sites + logical, optional :: sites_outside ! were there sites outside mesh? + + integer :: i ! loop index for local fission sites + integer :: ijk(3) ! indices on mesh + real(8) :: weight ! accumulated weight of sites + logical :: in_mesh ! was single site outside mesh? + logical :: outside ! was any site outside mesh? +#ifdef MPI + integer :: n ! total size of count variable +#endif + + ! initialize variables + count = ZERO + weight = ZERO + outside = .false. + + ! loop over fission sites and count how many are in each mesh box + FISSION_SITES: do i = 1, int(n_bank,4) + ! determine scoring bin for entropy mesh + call get_mesh_indices(m, fission_bank(i) % xyz, ijk, in_mesh) + + ! if outside mesh, skip particle + if (.not. in_mesh) then + outside = .true. + cycle + end if + + ! add weight + weight = weight + ONE + + ! add to appropriate mesh box + ! TODO: if tracking weight through bank, add weight instead + count(ijk(1),ijk(2),ijk(3)) = count(ijk(1),ijk(2),ijk(3)) + 1 + end do FISSION_SITES + +#ifdef MPI + ! determine total number of mesh cells + n = size(count,1) * size(count,2) * size(count,3) + + ! collect values from all processors + if (master) then + call MPI_REDUCE(MPI_IN_PLACE, count, n, MPI_REAL8, MPI_SUM, 0, & + MPI_COMM_WORLD, mpi_err) + else + call MPI_REDUCE(count, count, n, MPI_REAL8, MPI_SUM, 0, & + MPI_COMM_WORLD, mpi_err) + end if + + ! Check if there were sites outside the mesh for any processor + if (present(sites_outside)) then + call MPI_REDUCE(outside, sites_outside, 1, MPI_LOGICAL, MPI_LOR, 0, & + MPI_COMM_WORLD, mpi_err) + end if + + ! determine total weight of bank sites + call MPI_REDUCE(weight, total, 1, MPI_REAL8, MPI_SUM, 0, & + MPI_COMM_WORLD, mpi_err) +#else + total = weight + sites_outside = outside +#endif + + end subroutine count_fission_sites + end module mesh