From de24c7e8edb11e7cc8979331f3122b0a6abb17bf Mon Sep 17 00:00:00 2001
From: Paul Romano <romano7@gmail.com>
Date: Wed, 14 Dec 2011 00:00:33 -0600
Subject: [PATCH] Initial routine for calculating Shannon entropy.

---
 src/DEPENDENCIES   |  4 +++
 src/OBJECTS        |  1 +
 src/global.F90     |  7 ++++
 src/intercycle.F90 | 80 ++++++++++++++++++++++++++++++++++++++++++++++
 src/main.F90       |  4 +++
 5 files changed, 96 insertions(+)
 create mode 100644 src/intercycle.F90

diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES
index f6585ba4dc..7b80cc9a0b 100644
--- a/src/DEPENDENCIES
+++ b/src/DEPENDENCIES
@@ -95,6 +95,9 @@ input_xml.o: xml-fortran/templates/materials_t.o
 input_xml.o: xml-fortran/templates/settings_t.o
 input_xml.o: xml-fortran/templates/tallies_t.o
 
+intercycle.o: global.o
+intercycle.o: error.o
+
 interpolation.o: constants.o
 interpolation.o: endf_header.o
 interpolation.o: error.o
@@ -105,6 +108,7 @@ interpolation.o: string.o
 main.o: constants.o
 main.o: global.o
 main.o: initialize.o
+main.o: intercycle.o
 main.o: mpi_routines.o
 main.o: output.o
 main.o: particle_header.o
diff --git a/src/OBJECTS b/src/OBJECTS
index e8df9f0f7e..09f29ba6f8 100644
--- a/src/OBJECTS
+++ b/src/OBJECTS
@@ -15,6 +15,7 @@ geometry.o \
 geometry_header.o \
 global.o \
 initialize.o \
+intercycle.o \
 interpolation.o \
 input_xml.o \
 main.o \
diff --git a/src/global.F90 b/src/global.F90
index 814313c71e..4cfb6a6758 100644
--- a/src/global.F90
+++ b/src/global.F90
@@ -111,6 +111,13 @@ module global
   real(8) :: keff = ONE
   real(8) :: keff_std
 
+  ! Shannon entropy
+  logical :: entropy_on = .false.
+  real(8) :: entropy                ! value of shannon entropy
+  real(8) :: entropy_lower_left(3)  ! lower-left corner for entropy box
+  real(8) :: entropy_upper_right(3) ! upper-right corner for entropy box
+  real(8), allocatable :: entropy_p(:,:,:)
+
   ! ============================================================================
   ! PARALLEL PROCESSING VARIABLES
 
diff --git a/src/intercycle.F90 b/src/intercycle.F90
new file mode 100644
index 0000000000..9b8cf25a66
--- /dev/null
+++ b/src/intercycle.F90
@@ -0,0 +1,80 @@
+module intercycle
+
+  use global
+  use error, only: warning
+
+#ifdef MPI
+  use mpi
+#endif
+
+contains
+
+!===============================================================================
+! SHANNON_ENTROPY calculates the Shannon entropy of the fission source
+! distribution to assess source convergence
+!===============================================================================
+
+  subroutine shannon_entropy()
+
+    integer :: i_bank         ! 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, 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
+
+    ! 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
+       ! determine number of boxes in each direction
+       n = ceiling((n_particles/20)**(1.0/3.0))
+       n_box = n*n*n
+
+       ! determine width
+       width = (entropy_upper_right - entropy_lower_left)/n
+
+       ! allocate p
+       allocate(entropy_p(n,n,n))
+    end if
+
+    ! initialize p
+    entropy_p = ZERO
+
+    ! loop over fission sites and count how many are in each mesh box
+    FISSION_SITES: do i_bank = 1, n_bank
+       ! determine indices for entropy mesh box
+       i = (fission_bank(i_bank) % xyz(1) - entropy_lower_left(1))/n
+       j = (fission_bank(i_bank) % xyz(2) - entropy_lower_left(2))/n
+       k = (fission_bank(i_bank) % xyz(3) - entropy_lower_left(3))/n
+
+       ! if outside mesh, skip particle
+       if (i < 1 .or. i > n .or. j < 1 .or. &
+            j > n .or. k < 1 .or. k > n) then
+          message = "Fission source site outside of entropy box."
+          call warning()
+          cycle
+       end if
+
+       ! add to appropriate mesh box
+       entropy_p(i,j,k) = entropy_p(i,j,k) + 1
+    end do FISSION_SITES
+
+    ! 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)
+#endif
+
+    ! sum values to obtain shannon entropy
+    if (master) entropy = sum(entropy_p * log(entropy_p)/log(2.0))
+
+  end subroutine shannon_entropy
+
+end module intercycle
diff --git a/src/main.F90 b/src/main.F90
index 35479e9ea1..d7b11eb54e 100644
--- a/src/main.F90
+++ b/src/main.F90
@@ -3,6 +3,7 @@ program main
   use constants
   use global
   use initialize,      only: initialize_run
+  use intercycle,      only: shannon_entropy
   use mpi_routines,    only: synchronize_bank
   use output,          only: write_message, header, print_runtime
   use particle_header, only: Particle
@@ -117,6 +118,9 @@ contains
           call timer_stop(time_ic_tallies)
        end if
 
+       ! Calculate shannon entropy
+       if (entropy_on) call shannon_entropy()
+
        ! Distribute fission bank across processors evenly
        call synchronize_bank(i_cycle)