crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Allow k_generation and entropy to grow beyond n_max_batches

Browse source
This commit is contained in:
Paul Romano 2017-11-09 11:17:42 -06:00
parent 7d00c0b5eb
commit be28e13220
6 changed files with 40 additions and 31 deletions
Download patch file Download diff file Expand all files Collapse all files

28
src/eigenvalue.F90
View file

@ -301,8 +301,8 @@ contains
subroutine shannon_entropy()
integer :: ent_idx ! entropy index
integer :: i ! index for mesh elements
real(8) :: entropy_gen ! entropy at this generation
logical :: sites_outside ! were there sites outside entropy box?
associate (m => meshes(index_entropy_mesh))
@ -320,14 +320,16 @@ contains
! Normalize to total weight of bank sites
entropy_p = entropy_p / sum(entropy_p)
ent_idx = current_gen + gen_per_batch*(current_batch - 1)
entropy(ent_idx) = ZERO
entropy_gen = ZERO
do i = 1, size(entropy_p, 2)
if (entropy_p(1,i) > ZERO) then
entropy(ent_idx) = entropy(ent_idx) - &
entropy_gen = entropy_gen - &
entropy_p(1,i) * log(entropy_p(1,i))/log(TWO)
end if
end do
! Add value to vector
call entropy % push_back(entropy_gen)
end if
end associate
end subroutine shannon_entropy
@ -340,7 +342,7 @@ contains
subroutine calculate_generation_keff()
integer :: i ! overall generation
real(8) :: keff_reduced
#ifdef MPI
integer :: mpi_err ! MPI error code
#endif
@ -348,20 +350,18 @@ contains
! Get keff for this generation by subtracting off the starting value
keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH) - keff_generation
! Determine overall generation
i = overall_generation()
#ifdef MPI
! Combine values across all processors
call MPI_ALLREDUCE(keff_generation, k_generation(i), 1, MPI_REAL8, &
call MPI_ALLREDUCE(keff_generation, keff_reduced, 1, MPI_REAL8, &
MPI_SUM, mpi_intracomm, mpi_err)
#else
k_generation(i) = keff_generation
keff_reduced = keff_generation
#endif
! Normalize single batch estimate of k
! TODO: This should be normalized by total_weight, not by n_particles
k_generation(i) = k_generation(i) / n_particles
keff_reduced = keff_reduced / n_particles
call k_generation % push_back(keff_reduced)
end subroutine calculate_generation_keff
@ -385,12 +385,12 @@ contains
if (n <= 0) then
! For inactive generations, use current generation k as estimate for next
! generation
keff = k_generation(i)
keff = k_generation % data(i)
else
! Sample mean of keff
k_sum(1) = k_sum(1) + k_generation(i)
k_sum(2) = k_sum(2) + k_generation(i)**2
k_sum(1) = k_sum(1) + k_generation % data(i)
k_sum(2) = k_sum(2) + k_generation % data(i)**2
! Determine mean
keff = k_sum(1) / n

7
src/input_xml.F90
View file

@ -861,10 +861,9 @@ contains
call get_node_value(node_base, "generations_per_batch", gen_per_batch)
end if
! Allocate array for batch keff and entropy
allocate(k_generation(n_max_batches*gen_per_batch))
allocate(entropy(n_max_batches*gen_per_batch))
entropy = ZERO
! Preallocate space for keff and entropy by generation
call k_generation % reserve(n_max_batches*gen_per_batch)
call entropy % initialize(n_max_batches*gen_per_batch)
! Get the trigger information for keff
if (check_for_node(node_base, "keff_trigger")) then

8
src/output.F90
View file

@ -384,11 +384,11 @@ contains
! write out information about batch and generation
write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') &
trim(to_str(current_batch)) // "/" // trim(to_str(current_gen))
write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') k_generation(i)
write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') k_generation % data(i)
! write out entropy info
if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
entropy(i)
entropy % data(i)
if (n > 1) then
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') &
@ -418,11 +418,11 @@ contains
write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') &
trim(to_str(current_batch)) // "/" // trim(to_str(gen_per_batch))
write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') &
k_generation(i)
k_generation % data(i)
! write out entropy info
if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
entropy(i)
entropy % data(i)
! write out accumulated k-effective if after first active batch
if (n > 1) then

1
src/simulation.F90
View file

@ -78,6 +78,7 @@ contains
! Handle restart runs
if (restart_run .and. current_batch <= restart_batch) then
call replay_batch_history()
retval = 0
return
end if

12
src/simulation_header.F90
View file

@ -3,6 +3,7 @@ module simulation_header
use bank_header
use constants
use settings, only: gen_per_batch
use stl_vector, only: VectorReal
implicit none
@ -31,7 +32,7 @@ module simulation_header
integer(8) :: current_work ! index in source bank of current history simulated
! Temporary k-effective values
real(8), allocatable :: k_generation(:) ! single-generation estimates of k
type(VectorReal) :: k_generation ! single-generation estimates of k
real(8) :: keff = ONE ! average k over active batches
real(8) :: keff_std ! standard deviation of average k
real(8) :: k_col_abs = ZERO ! sum over batches of k_collision * k_absorption
@ -39,7 +40,7 @@ module simulation_header
real(8) :: k_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength
! Shannon entropy
real(8), allocatable :: entropy(:) ! shannon entropy at each generation
type(VectorReal) :: entropy ! shannon entropy at each generation
real(8), allocatable :: entropy_p(:,:) ! % of source sites in each cell
! Uniform fission source weighting
@ -85,11 +86,14 @@ contains
subroutine free_memory_simulation()
if (allocated(overlap_check_cnt)) deallocate(overlap_check_cnt)
if (allocated(k_generation)) deallocate(k_generation)
if (allocated(entropy)) deallocate(entropy)
if (allocated(entropy_p)) deallocate(entropy_p)
if (allocated(source_frac)) deallocate(source_frac)
if (allocated(work_index)) deallocate(work_index)
call k_generation % clear()
call k_generation % shrink_to_fit()
call entropy % clear()
call entropy % shrink_to_fit()
end subroutine free_memory_simulation
end module simulation_header

15
src/state_point.F90
View file

@ -121,8 +121,11 @@ contains
if (run_mode == MODE_EIGENVALUE) then
call write_dataset(file_id, "n_inactive", n_inactive)
call write_dataset(file_id, "generations_per_batch", gen_per_batch)
call write_dataset(file_id, "k_generation", k_generation)
call write_dataset(file_id, "entropy", entropy)
k = k_generation % size()
call write_dataset(file_id, "k_generation", k_generation % data(1:k))
if (entropy_on) then
call write_dataset(file_id, "entropy", entropy % data(1:k))
end if
call write_dataset(file_id, "k_col_abs", k_col_abs)
call write_dataset(file_id, "k_col_tra", k_col_tra)
call write_dataset(file_id, "k_abs_tra", k_abs_tra)
@ -707,10 +710,12 @@ contains
if (run_mode == MODE_EIGENVALUE) then
call read_dataset(int_array(1), file_id, "n_inactive")
call read_dataset(gen_per_batch, file_id, "generations_per_batch")
call read_dataset(k_generation(1:restart_batch*gen_per_batch), &
call read_dataset(k_generation % data(1:restart_batch*gen_per_batch), &
file_id, "k_generation")
call read_dataset(entropy(1:restart_batch*gen_per_batch), &
file_id, "entropy")
if (entropy_on) then
call read_dataset(entropy % data(1:restart_batch*gen_per_batch), &
file_id, "entropy")
end if
call read_dataset(k_col_abs, file_id, "k_col_abs")
call read_dataset(k_col_tra, file_id, "k_col_tra")
call read_dataset(k_abs_tra, file_id, "k_abs_tra")