From 1da3280eba23c73789858885393ffdfbef52d9a9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 14 Dec 2011 10:50:24 -0600 Subject: [PATCH] Finished Shannon entropy implementation. Closes github #17. --- src/DEPENDENCIES | 1 + src/intercycle.F90 | 131 ++++++++++++++++++++++++++++++++++++++++----- src/main.F90 | 12 +++-- src/tally.F90 | 79 +-------------------------- 4 files changed, 129 insertions(+), 94 deletions(-) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 7b80cc9a0b..8190c78b61 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -97,6 +97,7 @@ input_xml.o: xml-fortran/templates/tallies_t.o intercycle.o: global.o intercycle.o: error.o +intercycle.o: output.o interpolation.o: constants.o interpolation.o: endf_header.o diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 74b728768e..95dcd7003e 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -1,7 +1,10 @@ module intercycle + use ISO_FORTRAN_ENV + use global - use error, only: warning + use error, only: warning + use output, only: write_message #ifdef MPI use mpi @@ -16,10 +19,11 @@ contains subroutine shannon_entropy() - integer :: m ! index for bank sites integer :: i ! x-index for entropy mesh integer :: j ! y-index for entropy mesh integer :: k ! z-index for entropy mesh + integer :: m ! index for bank sites + integer(8) :: total_bank ! total # of fission bank sites integer, save :: n_box ! total # of boxes on mesh integer, save :: n ! # of boxes in each dimension real(8), save :: width(3) ! width of box in each dimension @@ -47,11 +51,11 @@ contains outside_box = .false. ! loop over fission sites and count how many are in each mesh box - FISSION_SITES: do m = 1, n_bank + FISSION_SITES: do m = 1, int(n_bank,4) ! determine indices for entropy mesh box - i = (fission_bank(m) % xyz(1) - entropy_lower_left(1))/width(1) + 1 - j = (fission_bank(m) % xyz(2) - entropy_lower_left(2))/width(2) + 1 - k = (fission_bank(m) % xyz(3) - entropy_lower_left(3))/width(3) + 1 + i = int((fission_bank(m) % xyz(1) - entropy_lower_left(1))/width(1)) + 1 + j = int((fission_bank(m) % xyz(2) - entropy_lower_left(2))/width(2)) + 1 + k = int((fission_bank(m) % xyz(3) - entropy_lower_left(3))/width(3)) + 1 ! if outside mesh, skip particle if (i < 1 .or. i > n .or. j < 1 .or. & @@ -70,18 +74,119 @@ contains call warning() end if - ! normalize to number of fission sites - entropy_p = entropy_p/n - - ! collect values from all processors #ifdef MPI - call MPI_REDUCE(MPI_IN_PLACE, entropy_p, n_box, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) + ! 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) entropy = sum(entropy_p * log(entropy_p)/log(2.0)) + if (master) then + entropy_p = entropy_p / total_bank + entropy = -sum(entropy_p * log(entropy_p)/log(2.0), entropy_p > ZERO) + end if end subroutine shannon_entropy +!=============================================================================== +! CALCULATE_KEFF calculates the single cycle estimate of keff as well as the +! mean and standard deviation of the mean for active cycles and displays them +!=============================================================================== + + subroutine calculate_keff(i_cycle) + + integer, intent(in) :: i_cycle ! index of current cycle + + integer(8) :: total_bank ! total number of source sites + integer :: n ! active cycle number + real(8) :: k_cycle ! single cycle estimate of keff + real(8), save :: k_sum ! accumulated keff + real(8), save :: k_sum_sq ! accumulated keff**2 + + message = "Calculate cycle keff..." + call write_message(8) + + ! initialize sum and square of sum at beginning of run + if (i_cycle == 1) then + k_sum = ZERO + k_sum_sq = ZERO + end if + +#ifdef MPI + ! Collect number bank sites onto master process + call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, & + & MPI_COMM_WORLD, mpi_err) +#else + total_bank = n_bank +#endif + + ! Collect statistics and print output + if (master) then + ! Since the creation of bank sites was originally weighted by the last + ! cycle keff, we need to multiply by that keff to get the current cycle's + ! value + + k_cycle = real(total_bank)/real(n_particles)*keff + + if (i_cycle > n_inactive) then + ! Active cycle number + n = i_cycle - n_inactive + + ! Accumulate cycle estimate of k + k_sum = k_sum + k_cycle + k_sum_sq = k_sum_sq + k_cycle*k_cycle + + ! Determine mean and standard deviation of mean + keff = k_sum/n + keff_std = sqrt((k_sum_sq/n - keff*keff)/n) + + ! Display output for this cycle + if (i_cycle > n_inactive+1) then + if (entropy_on) then + write(UNIT=OUTPUT_UNIT, FMT=103) i_cycle, k_cycle, entropy, & + keff, keff_std + else + write(UNIT=OUTPUT_UNIT, FMT=101) i_cycle, k_cycle, keff, keff_std + end if + else + if (entropy_on) then + write(UNIT=OUTPUT_UNIT, FMT=102) i_cycle, k_cycle, entropy + else + write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle + end if + end if + else + ! Display output for inactive cycle + if (entropy_on) then + write(UNIT=OUTPUT_UNIT, FMT=102) i_cycle, k_cycle, entropy + else + write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle + end if + keff = k_cycle + end if + end if + +#ifdef MPI + ! Broadcast new keff value to all processors + call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err) +#endif + +100 format (2X,I5,2X,F8.5) +101 format (2X,I5,2X,F8.5,9X,F8.5," +/-",F8.5) +102 format (2X,I5,2X,F8.5,3X,F8.5) +103 format (2X,I5,2X,F8.5,3X,F8.5,3X,F8.5," +/-",F8.5) + + end subroutine calculate_keff + end module intercycle diff --git a/src/main.F90 b/src/main.F90 index d7b11eb54e..c9534c09e1 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -3,7 +3,7 @@ program main use constants use global use initialize, only: initialize_run - use intercycle, only: shannon_entropy + use intercycle, only: shannon_entropy, calculate_keff use mpi_routines, only: synchronize_bank use output, only: write_message, header, print_runtime use particle_header, only: Particle @@ -13,7 +13,7 @@ program main use source, only: get_source_particle use string, only: int_to_str use tally, only: synchronize_tallies, write_tallies, & - tally_statistics, calculate_keff + tally_statistics use timing, only: timer_start, timer_stop #ifdef MPI @@ -64,6 +64,12 @@ contains tallies_on = .false. call timer_start(time_inactive) + ! Display column titles + message = " Cycle k(cycle) Entropy Average k" + call write_message(1) + message = " ===== ======== ======= ===================" + call write_message(1) + ! ========================================================================== ! LOOP OVER CYCLES CYCLE_LOOP: do i_cycle = 1, n_cycles @@ -73,7 +79,7 @@ contains message = "Simulating cycle " // trim(int_to_str(i_cycle)) // "..." call write_message(8) - + ! Set all tallies to zero n_bank = 0 diff --git a/src/tally.F90 b/src/tally.F90 index 29f144ac3e..8b039c8af6 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1,13 +1,11 @@ module tally - use ISO_FORTRAN_ENV - use constants use error, only: fatal_error use global use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices use mesh_header, only: StructuredMesh - use output, only: write_message, header + use output, only: header use search, only: binary_search use string, only: int_to_str, real_to_str use tally_header, only: TallyScore, TallyMapItem, TallyMapElement @@ -23,81 +21,6 @@ module tally contains -!=============================================================================== -! CALCULATE_KEFF calculates the single cycle estimate of keff as well as the -! mean and standard deviation of the mean for active cycles and displays them -!=============================================================================== - - subroutine calculate_keff(i_cycle) - - integer, intent(in) :: i_cycle ! index of current cycle - - integer(8) :: total_bank ! total number of source sites - integer :: n ! active cycle number - real(8) :: k_cycle ! single cycle estimate of keff - real(8), save :: k_sum ! accumulated keff - real(8), save :: k_sum_sq ! accumulated keff**2 - - message = "Calculate cycle keff..." - call write_message(8) - - ! initialize sum and square of sum at beginning of run - if (i_cycle == 1) then - k_sum = ZERO - k_sum_sq = ZERO - end if - -#ifdef MPI - ! Collect number bank sites onto master process - call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, & - & MPI_COMM_WORLD, mpi_err) -#else - total_bank = n_bank -#endif - - ! Collect statistics and print output - if (master) then - ! Since the creation of bank sites was originally weighted by the last - ! cycle keff, we need to multiply by that keff to get the current cycle's - ! value - - k_cycle = real(total_bank)/real(n_particles)*keff - - if (i_cycle > n_inactive) then - ! Active cycle number - n = i_cycle - n_inactive - - ! Accumulate cycle estimate of k - k_sum = k_sum + k_cycle - k_sum_sq = k_sum_sq + k_cycle*k_cycle - - ! Determine mean and standard deviation of mean - keff = k_sum/n - keff_std = sqrt((k_sum_sq/n - keff*keff)/n) - - ! Display output for this cycle - if (i_cycle > n_inactive+1) then - write(UNIT=OUTPUT_UNIT, FMT=101) i_cycle, k_cycle, keff, keff_std - else - write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle - end if - else - ! Display output for inactive cycle - write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle - keff = k_cycle - end if - end if - -#ifdef MPI - ! Broadcast new keff value to all processors - call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err) -#endif - -100 format (2X,I5,2X,F8.5) -101 format (2X,I5,2X,F8.5,9X,F8.5,1X,F8.5) - - end subroutine calculate_keff - !=============================================================================== ! CREATE_TALLY_MAP 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