Move code related to k_generation to C++

This commit is contained in:
Paul Romano 2018-10-12 08:01:23 -05:00
parent d7da87ca1d
commit 8ebd26b160
11 changed files with 126 additions and 102 deletions

View file

@ -408,10 +408,10 @@ constexpr int STANDARD_DEVIATION {3};
// Global tally parameters
constexpr int N_GLOBAL_TALLIES {4};
constexpr int K_COLLISION {1};
constexpr int K_ABSORPTION {2};
constexpr int K_TRACKLENGTH {3};
constexpr int LEAKAGE {4};
constexpr int K_COLLISION {0};
constexpr int K_ABSORPTION {1};
constexpr int K_TRACKLENGTH {2};
constexpr int LEAKAGE {3};
// Differential tally independent variables
constexpr int DIFF_DENSITY {1};

View file

@ -14,6 +14,7 @@ namespace openmc {
// Global variables
//==============================================================================
extern double keff_generation; //!< Single-generation k on each processor
extern std::vector<double> entropy; //!< Shannon entropy at each generation
extern xt::xtensor<double, 1> source_frac; //!< Source fraction for UFS
@ -24,6 +25,9 @@ extern "C" int64_t n_bank;
// Non-member functions
//==============================================================================
//! Collect/normalize the tracklength keff from each process
extern "C" void calculate_generation_keff();
//! Sample/redistribute source sites from accumulated fission sites
extern "C" void synchronize_bank();

View file

@ -56,6 +56,9 @@ void calculate_work();
//! Initialize simulation
extern "C" void openmc_simulation_init_c();
//! Initialize a fission generation
extern "C" void initialize_generation();
//! Determine overall generation number
extern "C" int overall_generation();

View file

@ -16,42 +16,15 @@ module eigenvalue
implicit none
real(8) :: keff_generation ! Single-generation k on each processor
real(8) :: k_sum(2) = ZERO ! Used to reduce sum and sum_sq
interface
subroutine calculate_generation_keff() bind(C)
end subroutine
end interface
contains
!===============================================================================
! CALCULATE_GENERATION_KEFF collects the single-processor tracklength k's onto
! the master processor and normalizes them. This should work whether or not the
! no-reduce method is being used.
!===============================================================================
subroutine calculate_generation_keff()
real(8) :: keff_reduced
#ifdef OPENMC_MPI
integer :: mpi_err ! MPI error code
#endif
! Get keff for this generation by subtracting off the starting value
keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH) - keff_generation
#ifdef OPENMC_MPI
! Combine values across all processors
call MPI_ALLREDUCE(keff_generation, keff_reduced, 1, MPI_REAL8, &
MPI_SUM, mpi_intracomm, mpi_err)
#else
keff_reduced = keff_generation
#endif
! Normalize single batch estimate of k
! TODO: This should be normalized by total_weight, not by n_particles
keff_reduced = keff_reduced / n_particles
call k_generation % push_back(keff_reduced)
end subroutine calculate_generation_keff
!===============================================================================
! CALCULATE_AVERAGE_KEFF calculates the mean and standard deviation of the mean
! of k-effective during active generations and broadcasts the mean to all
@ -76,12 +49,12 @@ contains
if (n <= 0) then
! For inactive generations, use current generation k as estimate for next
! generation
keff = k_generation % data(i)
keff = k_generation(i)
else
! Sample mean of keff
k_sum(1) = k_sum(1) + k_generation % data(i)
k_sum(2) = k_sum(2) + k_generation % data(i)**2
k_sum(1) = k_sum(1) + k_generation(i)
k_sum(2) = k_sum(2) + k_generation(i)**2
! Determine mean
keff = k_sum(1) / n

View file

@ -5,6 +5,7 @@
#include "xtensor/xview.hpp"
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/hdf5_interface.h"
#include "openmc/mesh.h"
@ -14,6 +15,7 @@
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/timer.h"
#include "openmc/tallies/tally.h"
#include <algorithm> // for min
#include <string>
@ -24,6 +26,7 @@ namespace openmc {
// Global variables
//==============================================================================
double keff_generation;
std::vector<double> entropy;
xt::xtensor<double, 1> source_frac;
@ -31,6 +34,28 @@ xt::xtensor<double, 1> source_frac;
// Non-member functions
//==============================================================================
void calculate_generation_keff()
{
auto gt = global_tallies();
// Get keff for this generation by subtracting off the starting value
keff_generation = gt(K_TRACKLENGTH, RESULT_VALUE) - keff_generation;
double keff_reduced;
#ifdef OPENMC_MPI
// Combine values across all processors
MPI_Allreduce(&keff_generation, &keff_reduced, 1, MPI_DOUBLE,
MPI_SUM, mpi::intracomm);
#else
keff_reduced = keff_generation;
#endif
// Normalize single batch estimate of k
// TODO: This should be normalized by total_weight, not by n_particles
keff_reduced /= settings::n_particles;
simulation::k_generation.push_back(keff_reduced);
}
void synchronize_bank()
{
time_bank.start();
@ -375,18 +400,34 @@ double ufs_get_weight(const Particle* p)
}
}
extern "C" void entropy_to_hdf5(hid_t group)
extern "C" void write_eigenvalue_hdf5(hid_t group)
{
write_dataset(group, "n_inactive", settings::n_inactive);
write_dataset(group, "generations_per_batch", settings::gen_per_batch);
write_dataset(group, "k_generation", simulation::k_generation);
if (settings::entropy_on) {
write_dataset(group, "entropy", entropy);
}
write_dataset(group, "k_col_abs", simulation::k_col_abs);
write_dataset(group, "k_col_tra", simulation::k_col_tra);
write_dataset(group, "k_abs_tra", simulation::k_abs_tra);
std::array<double, 2> k_combined;
openmc_get_keff(k_combined.data());
write_dataset(group, "k_combined", k_combined);
}
extern "C" void entropy_from_hdf5(hid_t group)
extern "C" void read_eigenvalue_hdf5(hid_t group)
{
read_dataset(group, "generations_per_batch", settings::gen_per_batch);
int n = simulation::restart_batch*settings::gen_per_batch;
simulation::k_generation.resize(n);
read_dataset(group, "k_generation", simulation::k_generation);
if (settings::entropy_on) {
read_dataset(group, "entropy", entropy);
}
read_dataset(group, "k_col_abs", simulation::k_col_abs);
read_dataset(group, "k_col_tra", simulation::k_col_tra);
read_dataset(group, "k_abs_tra", simulation::k_abs_tra);
}
extern "C" double entropy_c(int i)

View file

@ -223,7 +223,7 @@ contains
if (run_mode == MODE_EIGENVALUE) then
! Preallocate space for keff and entropy by generation
call k_generation % reserve(n_max_batches*gen_per_batch)
call k_generation_reserve(n_max_batches*gen_per_batch)
end if
! Particle tracks

View file

@ -339,7 +339,7 @@ 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 % data(i)
write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') k_generation(i)
! write out entropy info
if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
@ -373,7 +373,7 @@ 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 % data(i)
k_generation(i)
! write out entropy info
if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &

View file

@ -11,7 +11,7 @@ module simulation
use cmfd_header, only: cmfd_on
use constants, only: ZERO
use eigenvalue, only: calculate_average_keff, calculate_generation_keff, &
keff_generation, k_sum
k_sum
#ifdef _OPENMP
use eigenvalue, only: join_bank_from_threads
#endif
@ -57,6 +57,9 @@ module simulation
subroutine initialize_source() bind(C)
end subroutine
subroutine initialize_generation() bind(C)
end subroutine
function sample_external_source() result(site) bind(C)
import Bank
type(Bank) :: site
@ -222,30 +225,6 @@ contains
end subroutine initialize_batch
!===============================================================================
! INITIALIZE_GENERATION
!===============================================================================
subroutine initialize_generation()
interface
subroutine ufs_count_sites() bind(C)
end subroutine
end interface
if (run_mode == MODE_EIGENVALUE) then
! Reset number of fission bank sites
n_bank = 0
! Count source sites if using uniform fission source weighting
if (ufs) call ufs_count_sites()
! Store current value of tracklength k
keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH)
end if
end subroutine initialize_generation
!===============================================================================
! FINALIZE_GENERATION
!===============================================================================
@ -438,7 +417,7 @@ contains
! Reset global variables -- this is done before loading state point (as that
! will potentially populate k_generation and entropy)
current_batch = 0
call k_generation % clear()
call k_generation_clear()
call entropy_clear()
need_depletion_rx = .false.

View file

@ -1,6 +1,7 @@
#include "openmc/simulation.h"
#include "openmc/capi.h"
#include "openmc/eigenvalue.h"
#include "openmc/message_passing.h"
#include "openmc/settings.h"
#include "openmc/tallies/tally.h"
@ -62,7 +63,7 @@ int thread_id; //!< ID of a given thread
} // namespace simulation
//==============================================================================
// Functions
// Non-member functions
//==============================================================================
void openmc_simulation_init_c()
@ -71,6 +72,20 @@ void openmc_simulation_init_c()
calculate_work();
}
void initialize_generation()
{
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
// Reset number of fission bank sites
n_bank = 0;
// Count source sites if using uniform fission source weighting
if (settings::ufs_on) ufs_count_sites();
// Store current value of tracklength k
keff_generation = global_tallies()(K_TRACKLENGTH, RESULT_VALUE);
}
}
int overall_generation()
{
using namespace simulation;
@ -139,4 +154,13 @@ void broadcast_triggers()
}
#endif
//==============================================================================
// Fortran compatibility
//==============================================================================
extern "C" double k_generation(int i) { return simulation::k_generation.at(i - 1); }
extern "C" int k_generation_size() { return simulation::k_generation.size(); }
extern "C" void k_generation_clear() { simulation::k_generation.clear(); }
extern "C" void k_generation_reserve(int i) { simulation::k_generation.reserve(i); }
} // namespace openmc

View file

@ -39,7 +39,6 @@ module simulation_header
! K-EIGENVALUE SIMULATION VARIABLES
! Temporary k-effective values
type(VectorReal) :: k_generation ! single-generation estimates of k
real(C_DOUBLE), bind(C) :: keff ! average k over active batches
real(C_DOUBLE), bind(C) :: keff_std ! standard deviation of average k
real(C_DOUBLE), bind(C) :: k_col_abs ! sum over batches of k_collision * k_absorption
@ -71,6 +70,25 @@ module simulation_header
import C_INT
integer(C_INT) :: gen
end function overall_generation
function k_generation(i) result(k) bind(C)
import C_DOUBLE, C_INT
integer(C_INT), value :: i
real(C_DOUBLE) :: k
end function
function k_generation_size() result(sz) bind(C)
import C_INT
integer(C_INT) :: sz
end function
subroutine k_generation_clear() bind(C)
end subroutine
subroutine k_generation_reserve(i) bind(C)
import C_INT
integer(C_INT), value :: i
end subroutine
end interface
contains
@ -83,8 +101,7 @@ contains
if (allocated(work_index)) deallocate(work_index)
call k_generation % clear()
call k_generation % shrink_to_fit()
call k_generation_clear()
call entropy_clear()
end subroutine free_memory_simulation

View file

@ -63,15 +63,13 @@ contains
integer(C_INT) :: err
logical :: write_source_
integer :: i, j, k
integer :: i, j
integer :: i_xs
integer, allocatable :: id_array(:)
integer(HID_T) :: file_id
integer(HID_T) :: cmfd_group, tallies_group, tally_group, &
filters_group, filter_group, derivs_group, &
deriv_group, runtime_group
integer(C_INT) :: ignored_err
real(C_DOUBLE) :: k_combined(2)
character(MAX_WORD_LEN), allocatable :: str_array(:)
character(C_CHAR), pointer :: string(:)
character(len=:, kind=C_CHAR), allocatable :: filename_
@ -83,7 +81,7 @@ contains
import HID_T
integer(HID_T), value :: group
end subroutine
subroutine entropy_to_hdf5(group) bind(C)
subroutine write_eigenvalue_hdf5(group) bind(C)
import HID_T
integer(HID_T), value :: group
end subroutine
@ -173,16 +171,7 @@ contains
! Write out information for eigenvalue run
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)
k = k_generation % size()
call write_dataset(file_id, "k_generation", k_generation % data(1:k))
call entropy_to_hdf5(file_id)
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)
ignored_err = openmc_get_keff(k_combined)
call write_dataset(file_id, "k_combined", k_combined)
call write_eigenvalue_hdf5(file_id)
! Write out CMFD info
if (cmfd_on) then
@ -513,7 +502,9 @@ contains
character(MAX_WORD_LEN) :: word
interface
subroutine entropy_from_hdf5() bind(C)
subroutine read_eigenvalue_hdf5(group) bind(C)
import HID_T
integer(HID_T), value :: group
end subroutine
end interface
@ -592,16 +583,7 @@ contains
! Read information specific to eigenvalue run
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")
n = restart_batch*gen_per_batch
call k_generation % resize(n)
call read_dataset(k_generation % data(1:n), file_id, "k_generation")
call entropy_from_hdf5()
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")
call read_eigenvalue_hdf5(file_id)
! Take maximum of statepoint n_inactive and input n_inactive
n_inactive = max(n_inactive, int_array(1))
@ -634,13 +616,14 @@ contains
! of active cycle or inactive cycle
if (restart_batch > n_inactive) then
do i = n_inactive + 1, restart_batch
k_sum(1) = k_sum(1) + k_generation % data(i)
k_sum(2) = k_sum(2) + k_generation % data(i)**2
k_sum(1) = k_sum(1) + k_generation(i)
k_sum(2) = k_sum(2) + k_generation(i)**2
end do
n = gen_per_batch*n_realizations
keff = k_sum(1) / n
else
keff = k_generation % data(n)
n = k_generation_size()
keff = k_generation(n)
end if
current_batch = restart_batch