User input changes: origin is now lower_left for tallies and mesh. entropy_box is now entropy with a few sub-elements. Shannon entropy calc uses StructuredMesh type.

2026-07-28 14:15:42 -04:00 · 2012-02-08 15:35:50 -05:00 · 2012-02-08 15:35:50 -05:00 · bd633353bb
commit bd633353bb
parent 8fc482a2a4
11 changed files with 121 additions and 69 deletions
--- a/src/DEPENDENCIES
+++ b/src/DEPENDENCIES
@ -115,6 +115,8 @@ input_xml.o: xml-fortran/templates/tallies_t.o

 intercycle.o: error.o
 intercycle.o: global.o
+intercycle.o: mesh.o
+intercycle.o: mesh_header.o
 intercycle.o: output.o
 intercycle.o: particle_header.o
 intercycle.o: random_lcg.o
--- a/src/geometry_header.F90
+++ b/src/geometry_header.F90
@ -27,8 +27,8 @@ module geometry_header
     integer :: level    ! Level of lattice
     integer :: n_x      ! number of lattice cells in x-direction
     integer :: n_y      ! number of lattice cells in y-direction
-     real(8) :: x0       ! x-coordinate of lattice origin
-     real(8) :: y0       ! y-coordinate of lattice origin
+     real(8) :: x0       ! x-coordinate of lower-left corner of lattice
+     real(8) :: y0       ! y-coordinate of lower-left corner of lattice
     real(8) :: width_x  ! width of lattice cell 
     real(8) :: width_y  ! width of lattice cell
     integer, allocatable :: element(:,:) ! specified universes
--- a/src/global.F90
+++ b/src/global.F90
@ -125,10 +125,9 @@ module global

  ! 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(:,:,:)
+  real(8) :: entropy                   ! value of shannon entropy
+  real(8), allocatable :: entropy_p(:) ! fraction of source sites in each cell
+  type(StructuredMesh), pointer :: entropy_mesh

  ! ============================================================================
  ! PARALLEL PROCESSING VARIABLES
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@ -145,24 +145,50 @@ contains
       trace_particle = trace_(2)
    end if

-    ! Entropy box
-    if (associated(entropy_box_)) then
+    ! Shannon Entropy mesh
+    if (size(entropy_) > 0) then
       ! Check to make sure enough values were supplied
-       if (size(entropy_box_) /= 6) then
-          message = "Need to supply lower-left and upper-right coordinates &
-               &for Shannon entropy box."
+       if (size(entropy_(1) % lower_left) /= 3) then
+          message = "Need to specify (x,y,z) coordinates of lower-left corner &
+               &of Shannon entropy mesh."
+       elseif (size(entropy_(1) % upper_right) /= 3) then
+          message = "Need to specify (x,y,z) coordinates of upper-right 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))
+
+       ! Copy values
+       entropy_mesh % lower_left  = entropy_(1) % lower_left
+       entropy_mesh % upper_right = entropy_(1) % upper_right
+
+       ! Check on values provided
+       if (.not. all(entropy_mesh % upper_right > entropy_mesh % lower_left)) then
+          message = "Upper-right coordinate must be greater than lower-left &
+               &coordinate for Shannon entropy mesh."
          call fatal_error()
       end if

-       ! Copy values
-       entropy_lower_left = entropy_box_(1:3)
-       entropy_upper_right = entropy_box_(4:6)
+       ! Check if dimensions were specified -- if not, they will be calculated
+       ! automatically upon first entry into shannon_entropy

-       ! Check on values provided
-       if (.not. all(entropy_upper_right > entropy_lower_left)) then
-          message = "Upper-right coordinate must be greater than lower-left &
-               &coordinate for Shannon entropy box."
-          call fatal_error()
+       if (associated(entropy_(1) % dimension)) then
+          ! If so, make sure proper number of values were given
+          if (size(entropy_(1) % dimension) /= 3) then
+             message = "Dimension of entropy mesh must be given as three &
+                  &integers."
+             call fatal_error()
+          end if
+
+          ! Allocate dimensions
+          entropy_mesh % n_dimension = 3
+          allocate(entropy_mesh % dimension(3))
+
+          ! Copy dimensions
+          entropy_mesh % dimension = entropy_(1) % dimension
       end if

       ! Turn on Shannon entropy calculation
@ -409,19 +435,19 @@ contains
       l % n_x = n_x
       l % n_y = n_y

-       ! Read lattice origin location
-       if (size(lattice_(i) % dimension) /= size(lattice_(i) % origin)) then
-          message = "Number of entries on <origin> must be the same as " // &
+       ! Read lattice lower-left location
+       if (size(lattice_(i) % dimension) /= size(lattice_(i) % lower_left)) then
+          message = "Number of entries on <lower_left> must be the same as " // &
               "the number of entries on <dimension>."
          call fatal_error()
       end if
-       l % x0 = lattice_(i) % origin(1)
-       l % y0 = lattice_(i) % origin(2)
+       l % x0 = lattice_(i) % lower_left(1)
+       l % y0 = lattice_(i) % lower_left(2)

       ! Read lattice widths
-       if (size(lattice_(i) % width) /= size(lattice_(i) % origin)) then
+       if (size(lattice_(i) % width) /= size(lattice_(i) % lower_left)) then
          message = "Number of entries on <width> must be the same as " // &
-               "the number of entries on <origin>."
+               "the number of entries on <lower_left>."
          call fatal_error()
       end if
       l % width_x = lattice_(i) % width(1)
@ -697,31 +723,31 @@ contains

       ! Allocate attribute arrays
       allocate(m % dimension(n))
-       allocate(m % origin(n))
+       allocate(m % lower_left(n))
       allocate(m % width(n))
       allocate(m % upper_right(n))

       ! Read dimensions in each direction
       m % dimension = mesh_(i) % dimension

-       ! Read mesh origin location
-       if (m % n_dimension /= size(mesh_(i) % origin)) then
-          message = "Number of entries on <origin> must be the same as the &
-               &number of entries on <dimension>."
+       ! Read mesh lower-left corner location
+       if (m % n_dimension /= size(mesh_(i) % lower_left)) then
+          message = "Number of entries on <lower_left> must be the same as &
+               &the number of entries on <dimension>."
          call fatal_error()
       end if
-       m % origin = mesh_(i) % origin
+       m % lower_left = mesh_(i) % lower_left

       ! Read mesh widths
-       if (size(mesh_(i) % width) /= size(mesh_(i) % origin)) then
+       if (size(mesh_(i) % width) /= size(mesh_(i) % lower_left)) then
          message = "Number of entries on <width> must be the same as the &
-               &number of entries on <origin>."
+               &number of entries on <lower_left>."
          call fatal_error()
       end if
       m % width = mesh_(i) % width

       ! Set upper right coordinate
-       m % upper_right = m % origin + m % dimension * m % width
+       m % upper_right = m % lower_left + m % dimension * m % width

       ! Add mesh to dictionary
       call dict_add_key(mesh_dict, m % id, i)
@ -885,8 +911,8 @@ contains
       t % filters = filters(1:n_filters)

       ! Read macro reactions
-       if (len_trim(tally_(i) % macros) > 0) then
-          call split_string(tally_(i) % macros, words, n_words)
+       if (len_trim(tally_(i) % scores) > 0) then
+          call split_string(tally_(i) % scores, words, n_words)
          allocate(t % score_bins(n_words))
          do j = 1, n_words
             word = words(j)
--- a/src/intercycle.F90
+++ b/src/intercycle.F90
@ -4,6 +4,8 @@ module intercycle

  use error,           only: fatal_error, warning
  use global
+  use mesh,            only: get_mesh_bin
+  use mesh_header,     only: StructuredMesh
  use output,          only: write_message
  use particle_header, only: Particle, initialize_particle
  use random_lcg,      only: prn, set_particle_seed, prn_skip
@ -328,15 +330,16 @@ contains

  subroutine shannon_entropy()

-    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 :: i              ! 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
-    integer, save :: n        ! # of boxes in each dimension
-    real(8), save :: width(3) ! width of box in each dimension
    logical :: outside_box    ! were there sites outside entropy box?
+    type(StructuredMesh), pointer :: m => null()
+
+    ! 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
@ -344,15 +347,27 @@ contains
    ! 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
+       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)**(1.0/3.0))

-       ! determine width
-       width = (entropy_upper_right - entropy_lower_left)/n
+          ! copy dimensions
+          m % n_dimension = 3
+          allocate(m % dimension(3))
+          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,n,n))
+       allocate(entropy_p(n_box))
    end if

    ! initialize p
@ -360,21 +375,18 @@ contains
    outside_box = .false.

    ! loop over fission sites and count how many are in each mesh box
-    FISSION_SITES: do m = 1, int(n_bank,4)
-       ! determine indices for entropy mesh box
-       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
+    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 (i < 1 .or. i > n .or. j < 1 .or. &
-            j > n .or. k < 1 .or. k > n) then
+       if (bin == NO_BIN_FOUND) then
          outside_box = .true.
          cycle
       end if

       ! add to appropriate mesh box
-       entropy_p(i,j,k) = entropy_p(i,j,k) + 1
+       entropy_p(bin) = entropy_p(bin) + 1
    end do FISSION_SITES

    ! display warning message if there were sites outside entropy box
@ -402,8 +414,15 @@ contains

    ! sum values to obtain shannon entropy
    if (master) then
+       ! Normalize to number of bank sites
       entropy_p = entropy_p / total_bank
-       entropy = -sum(entropy_p * log(entropy_p)/log(2.0), entropy_p > ZERO)
+
+       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
+       end do
    end if

  end subroutine shannon_entropy
--- a/src/mesh.F90
+++ b/src/mesh.F90
@ -26,13 +26,13 @@ contains
    n = m % n_dimension

    ! Check for cases where particle is outside of mesh
-    if (xyz(1) < m % origin(1)) then
+    if (xyz(1) < m % lower_left(1)) then
       bin = NO_BIN_FOUND
       return
    elseif (xyz(1) > m % upper_right(1)) then
       bin = NO_BIN_FOUND
       return
-    elseif (xyz(2) < m % origin(2)) then
+    elseif (xyz(2) < m % lower_left(2)) then
       bin = NO_BIN_FOUND
       return
    elseif (xyz(2) > m % upper_right(2)) then
@ -40,7 +40,7 @@ contains
       return
    end if
    if (n > 2) then
-       if (xyz(3) < m % origin(3)) then
+       if (xyz(3) < m % lower_left(3)) then
          bin = NO_BIN_FOUND
          return
       elseif (xyz(3) > m % upper_right(3)) then
@ -73,7 +73,7 @@ contains
    logical, intent(out)          :: in_mesh

    ! Find particle in mesh
-    ijk = ceiling((xyz - m % origin)/m % width)
+    ijk = ceiling((xyz - m % lower_left)/m % width)

    ! Determine if particle is in mesh
    if (any(ijk < 1) .or. any(ijk > m % dimension)) then
--- a/src/mesh_header.F90
+++ b/src/mesh_header.F90
@ -12,7 +12,7 @@ module mesh_header
     integer :: type
     integer :: n_dimension
     integer, allocatable :: dimension(:)
-     real(8), allocatable :: origin(:)
+     real(8), allocatable :: lower_left(:)
     real(8), allocatable :: upper_right(:)
     real(8), allocatable :: width(:)
  end type StructuredMesh
--- a/src/tally.F90
+++ b/src/tally.F90
@ -883,9 +883,9 @@ contains
       ! Bounding coordinates
       do j = 1, 3
          if (uvw(j) > 0) then
-             xyz_cross(j) = m % origin(j) + ijk0(j) * m % width(j)
+             xyz_cross(j) = m % lower_left(j) + ijk0(j) * m % width(j)
          else
-             xyz_cross(j) = m % origin(j) + (ijk0(j) - 1) * m % width(j)
+             xyz_cross(j) = m % lower_left(j) + (ijk0(j) - 1) * m % width(j)
          end if
       end do

--- a/src/xml-fortran/templates/geometry_t.xml
+++ b/src/xml-fortran/templates/geometry_t.xml
@ -24,7 +24,7 @@
    <component name="id" type="integer" />
    <component name="type" type="word" length="12" default="'rectangular'" />
    <component name="dimension" type="integer-array" />
-    <component name="origin" type="double-array" />
+    <component name="lower_left" type="double-array" />
    <component name="width" type="double-array" />
    <component name="universes" type="integer-array" />
  </typedef>
--- a/src/xml-fortran/templates/settings_t.xml
+++ b/src/xml-fortran/templates/settings_t.xml
@ -19,6 +19,12 @@
    <component name="weight_avg" type="double" default="1.0" />
  </typedef>

+  <typedef name="entropy_xml">
+    <component name="dimension" type="integer-array" />
+    <component name="lower_left" type="double-array" />
+    <component name="upper_right" type="double-array" />
+  </typedef>
+
  <variable name="criticality" type="criticality_xml" />
  <variable name="verbosity_" tag="verbosity" type="integer" />
  <variable name="source_" tag="source" type="source_xml" />
@ -27,6 +33,6 @@
  <variable name="cutoff_" tag="cutoff" type="cutoff_xml" dimension="1" />
  <variable name="cross_sections_" tag="cross_sections" type="word" length="255" />
  <variable name="trace_" tag="trace" type="integer-array" />
-  <variable name="entropy_box_" tag="entropy_box" type="double-array" />
+  <variable name="entropy_" tag="entropy" type="entropy_xml" dimension="1" />

 </template>
--- a/src/xml-fortran/templates/tallies_t.xml
+++ b/src/xml-fortran/templates/tallies_t.xml
@ -7,7 +7,7 @@
    <component name="id" type="integer" />
    <component name="type" type="word" length="12" />
    <component name="dimension" type="integer-array" />
-    <component name="origin" type="double-array" />
+    <component name="lower_left" type="double-array" />
    <component name="width" type="double-array" />
  </typedef>

@ -25,7 +25,7 @@
  <typedef name="tally_xml">
    <component name="id" type="integer" />
    <component name="filters" type="filter_xml" />
-    <component name="macros" type="word" length="250" default="''" />
+    <component name="scores" type="word" length="250" default="''" />
    <component name="reactions" type="word" length="250" default="''" />
    <component name="nuclides" type="word" length="250" default="''" />
  </typedef>