mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
Finished Shannon entropy implementation. Closes github #17.
This commit is contained in:
Paul Romano
parent
69f8f9bd46
commit
1da3280eba
4 changed files with 129 additions and 94 deletions
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
12
src/main.F90
12
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
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue