From be28e1322031e04969f26f363d021ba4411de9d8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 9 Nov 2017 11:17:42 -0600 Subject: [PATCH] Allow k_generation and entropy to grow beyond n_max_batches --- src/eigenvalue.F90 | 28 ++++++++++++++-------------- src/input_xml.F90 | 7 +++---- src/output.F90 | 8 ++++---- src/simulation.F90 | 1 + src/simulation_header.F90 | 12 ++++++++---- src/state_point.F90 | 15 ++++++++++----- 6 files changed, 40 insertions(+), 31 deletions(-) diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 5f99d1580a..27ae4063ea 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 71ebe3ae6b..69b519ea5a 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 diff --git a/src/output.F90 b/src/output.F90 index f820ce4b9d..3d02489f78 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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 diff --git a/src/simulation.F90 b/src/simulation.F90 index ad2841f54d..e3eba0d512 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -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 diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index cfb75b4492..9719cc72e3 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -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 diff --git a/src/state_point.F90 b/src/state_point.F90 index 265cc49c02..e93d6e4729 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -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")