From bd633353bb86046a2bc8ff848ae80decf3c9f46b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 8 Feb 2012 15:35:50 -0500 Subject: [PATCH] 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. --- src/DEPENDENCIES | 2 + src/geometry_header.F90 | 4 +- src/global.F90 | 7 +- src/input_xml.F90 | 88 +++++++++++++++--------- src/intercycle.F90 | 61 ++++++++++------ src/mesh.F90 | 8 +-- src/mesh_header.F90 | 2 +- src/tally.F90 | 4 +- src/xml-fortran/templates/geometry_t.xml | 2 +- src/xml-fortran/templates/settings_t.xml | 8 ++- src/xml-fortran/templates/tallies_t.xml | 4 +- 11 files changed, 121 insertions(+), 69 deletions(-) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index e5d327bb8a..8310c8ab04 100644 --- 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 diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 35de1996bf..fbcf0353a9 100644 --- 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 diff --git a/src/global.F90 b/src/global.F90 index 3c6fd8c1b2..5c09afe46f 100644 --- 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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 7fe7a1b4a2..ad2f3d60ac 100644 --- 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 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 must be the same as " // & "the number of entries on ." 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 must be the same as " // & - "the number of entries on ." + "the number of entries on ." 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 must be the same as the & - &number of entries on ." + ! Read mesh lower-left corner location + if (m % n_dimension /= size(mesh_(i) % lower_left)) then + message = "Number of entries on must be the same as & + &the number of entries on ." 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 must be the same as the & - &number of entries on ." + &number of entries on ." 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) diff --git a/src/intercycle.F90 b/src/intercycle.F90 index e95cda1419..bfe1e71cc1 100644 --- 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 diff --git a/src/mesh.F90 b/src/mesh.F90 index 8e2de583b4..69f7b47eaf 100644 --- 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 diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 index ff4d30cf37..b70bbd45b4 100644 --- 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 diff --git a/src/tally.F90 b/src/tally.F90 index bab721bf66..eda14464e1 100644 --- 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 diff --git a/src/xml-fortran/templates/geometry_t.xml b/src/xml-fortran/templates/geometry_t.xml index 6e0a7a368a..ef56586523 100644 --- a/src/xml-fortran/templates/geometry_t.xml +++ b/src/xml-fortran/templates/geometry_t.xml @@ -24,7 +24,7 @@ - + diff --git a/src/xml-fortran/templates/settings_t.xml b/src/xml-fortran/templates/settings_t.xml index 0bfa90278c..3f5456da29 100644 --- a/src/xml-fortran/templates/settings_t.xml +++ b/src/xml-fortran/templates/settings_t.xml @@ -19,6 +19,12 @@ + + + + + + @@ -27,6 +33,6 @@ - + diff --git a/src/xml-fortran/templates/tallies_t.xml b/src/xml-fortran/templates/tallies_t.xml index d26545606d..d88937f5a9 100644 --- a/src/xml-fortran/templates/tallies_t.xml +++ b/src/xml-fortran/templates/tallies_t.xml @@ -7,7 +7,7 @@ - + @@ -25,7 +25,7 @@ - +