crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge remote-tracking branch 'upstream/develop' into cpp_tallies

Browse source
This commit is contained in:
Sterling Harper 2018-10-28 22:28:18 -04:00
commit f7c45a3fd1
47 changed files with 1220 additions and 1308 deletions
Download patch file Download diff file Expand all files Collapse all files

3
include/openmc/capi.h
View file

@ -95,7 +95,7 @@ extern "C" {
int openmc_sphharm_filter_get_cosine(int32_t index, char cosine[]);
int openmc_sphharm_filter_set_order(int32_t index, int order);
int openmc_sphharm_filter_set_cosine(int32_t index, const char cosine[]);
int openmc_statepoint_write(const char filename[]);
int openmc_statepoint_write(const char filename[], bool* write_source);
int openmc_tally_allocate(int32_t index, const char* type);
int openmc_tally_get_active(int32_t index, bool* active);
int openmc_tally_get_estimator(int32_t index, int32_t* estimator);
@ -136,7 +136,6 @@ extern "C" {
// Global variables
extern char openmc_err_msg[256];
extern int32_t n_cells;
extern int32_t n_filters;
extern int32_t n_lattices;
extern int32_t n_materials;
extern int n_nuclides;

17
include/openmc/container_util.h Normal file
View file

@ -0,0 +1,17 @@
#ifndef OPENMC_CONTAINER_UTIL_H
#define OPENMC_CONTAINER_UTIL_H
#include <algorithm> // for find
#include <iterator> // for begin, end
namespace openmc {
template<class C, class T>
inline bool contains(const C& v, const T& x)
{
return std::end(v) != std::find(std::begin(v), std::end(v), x);
}
}
#endif // OPENMC_CONTAINER_UTIL_H

27
include/openmc/eigenvalue.h
View file

@ -1,6 +1,10 @@
//! \file eigenvalue.h
//! \brief Data/functions related to k-eigenvalue calculations
#ifndef OPENMC_EIGENVALUE_H
#define OPENMC_EIGENVALUE_H
#include <array>
#include <cstdint> // for int64_t
#include <vector>
@ -15,6 +19,7 @@ namespace openmc {
//==============================================================================
extern double keff_generation; //!< Single-generation k on each processor
extern std::array<double, 2> k_sum; //!< Used to reduce sum and sum_sq
extern std::vector<double> entropy; //!< Shannon entropy at each generation
extern xt::xtensor<double, 1> source_frac; //!< Source fraction for UFS
@ -28,6 +33,28 @@ extern "C" int64_t n_bank;
//! Collect/normalize the tracklength keff from each process
extern "C" void calculate_generation_keff();
//! Calculate mean/standard deviation of keff during active generations
//!
//! This function sets the global variables keff and keff_std which represent
//! the mean and standard deviation of the mean of k-effective over active
//! generations. It also broadcasts the value from the master process.
extern "C" void calculate_average_keff();
//! Calculates a minimum variance estimate of k-effective
//!
//! The minimum variance estimate is based on a linear combination of the
//! collision, absorption, and tracklength estimates. The theory behind this can
//! be found in M. Halperin, "Almost linearly-optimum combination of unbiased
//! estimates," J. Am. Stat. Assoc., 56, 36-43 (1961),
//! doi:10.1080/01621459.1961.10482088. The implementation here follows that
//! described in T. Urbatsch et al., "Estimation and interpretation of keff
//! confidence intervals in MCNP," Nucl. Technol., 111, 169-182 (1995).
//!
//! \param[out] k_combined Estimate of k-effective and its standard deviation
//! \return Error status
extern "C" int openmc_get_keff(double* k_combined);
//! Sample/redistribute source sites from accumulated fission sites
extern "C" void synchronize_bank();

2
include/openmc/finalize.h
View file

@ -3,8 +3,6 @@
namespace openmc {
extern "C" void openmc_free_bank();
} // namespace openmc
#endif // OPENMC_FINALIZE_H

10
include/openmc/hdf5_interface.h
View file

@ -52,7 +52,7 @@ void write_string(hid_t group_id, const char* name, const std::string& buffer,
std::vector<hsize_t> attribute_shape(hid_t obj_id, const char* name);
std::vector<std::string> dataset_names(hid_t group_id);
void ensure_exists(hid_t group_id, const char* name);
void ensure_exists(hid_t obj_id, const char* name, bool attribute=false);
std::vector<std::string> group_names(hid_t group_id);
std::vector<hsize_t> object_shape(hid_t obj_id);
std::string object_name(hid_t obj_id);
@ -202,7 +202,13 @@ read_attribute(hid_t obj_id, const char* name, std::vector<std::string>& vec)
read_attr_string(obj_id, name, n, buffer[0]);
for (int i = 0; i < m; ++i) {
vec.emplace_back(&buffer[i][0], std::min(strlen(buffer[i]), n));
// Determine proper length of string -- strlen doesn't work because
// buffer[i] might not have any null characters
std::size_t k = 0;
for (; k < n; ++k) if (buffer[i][k] == '\0') break;
// Create string based on (char*, size_t) constructor
vec.emplace_back(&buffer[i][0], k);
}
}

3
include/openmc/initialize.h
View file

@ -5,9 +5,6 @@
#include "mpi.h"
#endif
extern "C" void print_usage();
extern "C" void print_version();
namespace openmc {
int parse_command_line(int argc, char* argv[]);

9
include/openmc/mgxs_interface.h
View file

@ -7,6 +7,7 @@
#include "hdf5_interface.h"
#include "mgxs.h"
#include <vector>
namespace openmc {
@ -17,9 +18,9 @@ namespace openmc {
extern std::vector<Mgxs> nuclides_MG;
extern std::vector<Mgxs> macro_xs;
extern "C" int num_energy_groups;
//TODO: When more of the Fortran is converted (input_xml, tallies, etc, also
// bring over energy_bin_avg, energy_bins, etc, as vectors)
extern std::vector<double> energy_bins;
extern std::vector<double> energy_bin_avg;
extern std::vector<double> rev_energy_bins;
//==============================================================================
// Mgxs data loading interface methods
@ -39,6 +40,8 @@ create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[],
int n_temps, const double temps[], const double atom_densities[],
double tolerance, int& method);
extern "C" void read_mg_cross_sections_header_c(hid_t file_id);
//==============================================================================
// Mgxs tracking/transport/tallying interface methods
//==============================================================================

2
include/openmc/output.h
View file

@ -20,7 +20,7 @@ void header(const char* msg, int level);
//! Display information regarding cell overlap checking.
//==============================================================================
extern "C" void print_overlap_check();
void print_overlap_check();
extern "C" void title();

7
include/openmc/particle.h
View file

@ -146,9 +146,7 @@ extern "C" {
//! site may have been produced from an external source, from fission, or
//! simply as a secondary particle.
//! \param src Source site data
//! \param run_CE Whether continuous-energy data is being used
//! \param energy_bin_avg An array of energy group bin averages
void from_source(const Bank* src, bool run_CE, const double* energy_bin_avg);
void from_source(const Bank* src);
//! mark a particle as lost and create a particle restart file
//! \param message A warning message to display
@ -174,8 +172,7 @@ extern "C" {
void particle_create_secondary(Particle* p, const double* uvw, double E,
int type, bool run_CE);
void particle_initialize(Particle* p);
void particle_from_source(Particle* p, const Bank* src, bool run_CE,
const double* energy_bin_avg);
void particle_from_source(Particle* p, const Bank* src);
void particle_mark_as_lost(Particle* p, const char* message);
void particle_write_restart(Particle* p);

1
include/openmc/random_lcg.h
View file

@ -17,6 +17,7 @@ extern "C" const int STREAM_SOURCE;
extern "C" const int STREAM_URR_PTABLE;
extern "C" const int STREAM_VOLUME;
extern "C" const int STREAM_PHOTON;
constexpr int64_t DEFAULT_SEED = 1;
//==============================================================================
//! Generate a pseudo-random number using a linear congruential generator.

10
include/openmc/settings.h
View file

@ -7,6 +7,8 @@
#include <array>
#include <cstdint>
#include <string>
#include <unordered_set>
#include <vector>
#include "pugixml.hpp"
@ -66,23 +68,25 @@ extern "C" int32_t gen_per_batch; //!< number of generations per batch
extern "C" int64_t n_particles; //!< number of particles per generation
extern "C" int electron_treatment; //!< how to treat secondary electrons
extern "C" double energy_cutoff[4]; //!< Energy cutoff in [eV] for each particle type
extern "C" std::array<double, 4> energy_cutoff; //!< Energy cutoff in [eV] for each particle type
extern "C" int legendre_to_tabular_points; //!< number of points to convert Legendres
extern "C" int max_order; //!< Maximum Legendre order for multigroup data
extern "C" int n_log_bins; //!< number of bins for logarithmic energy grid
extern "C" int n_max_batches; //!< Maximum number of batches
extern "C" int res_scat_method; //!< resonance upscattering method
extern "C" double res_scat_energy_min; //!< Min energy in [eV] for res. upscattering
extern "C" double res_scat_energy_max; //!< Max energy in [eV] for res. upscattering
extern "C" int run_mode; //!< Run mode (eigenvalue, fixed src, etc.)
extern std::unordered_set<int> sourcepoint_batch; //!< Batches when source should be written
extern std::unordered_set<int> statepoint_batch; //!< Batches when state should be written
extern "C" int temperature_method; //!< method for choosing temperatures
extern "C" double temperature_tolerance; //!< Tolerance in [K] on choosing temperatures
extern "C" double temperature_default; //!< Default T in [K]
extern "C" double temperature_range[2]; //!< Min/max T in [K] over which to load xs
extern "C" std::array<double, 2> temperature_range; //!< Min/max T in [K] over which to load xs
extern "C" int trace_batch; //!< Batch to trace particle on
extern "C" int trace_gen; //!< Generation to trace particle on
extern "C" int64_t trace_particle; //!< Particle ID to enable trace on
extern std::vector<std::array<int, 3>> track_identifiers; //!< Particle numbers for writing tracks
extern "C" int trigger_batch_interval; //!< Batch interval for triggers
extern "C" int verbosity; //!< How verbose to make output
extern "C" double weight_cutoff; //!< Weight cutoff for Russian roulette

29
include/openmc/simulation.h
View file

@ -4,11 +4,17 @@
#ifndef OPENMC_SIMULATION_H
#define OPENMC_SIMULATION_H
#include "openmc/particle.h"
#include <cstdint>
#include <vector>
namespace openmc {
constexpr int STATUS_EXIT_NORMAL {0};
constexpr int STATUS_EXIT_MAX_BATCH {1};
constexpr int STATUS_EXIT_ON_TRIGGER {2};
//==============================================================================
// Global variable declarations
//==============================================================================
@ -18,6 +24,7 @@ namespace simulation {
extern "C" int current_batch; //!< current batch
extern "C" int current_gen; //!< current fission generation
extern "C" int64_t current_work; //!< index in source back of current particle
extern "C" bool initialized; //!< has simulation been initialized?
extern "C" double keff; //!< average k over batches
extern "C" double keff_std; //!< standard deviation of average k
extern "C" double k_col_abs; //!< sum over batches of k_collision * k_absorption
@ -28,7 +35,6 @@ extern "C" int n_lost_particles; //!< cumulative number of lost particles
extern "C" bool need_depletion_rx; //!< need to calculate depletion rx?
extern "C" int restart_batch; //!< batch at which a restart job resumed
extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied?
extern "C" bool simulation_initialized; //!< has simulation been initialized?
extern "C" int total_gen; //!< total number of generations simulated
extern "C" int64_t work; //!< number of particles per process
@ -53,18 +59,29 @@ extern "C" int thread_id; //!< ID of a given thread
//! Determine number of particles to transport per process
void calculate_work();
//! Initialize simulation
extern "C" void openmc_simulation_init_c();
//! Initialize a batch
void initialize_batch();
//! Initialize a fission generation
extern "C" void initialize_generation();
void initialize_generation();
void initialize_history(Particle* p, int64_t index_source);
//! Finalize a batch
//!
//! Handles synchronization and accumulation of tallies, calculation of Shannon
//! entropy, getting single-batch estimate of keff, and turning on tallies when
//! appropriate
void finalize_batch();
//! Finalize a fission generation
void finalize_generation();
//! Determine overall generation number
extern "C" int overall_generation();
#ifdef OPENMC_MPI
extern "C" void broadcast_results();
extern "C" void broadcast_triggers();
void broadcast_results();
#endif
} // namespace openmc

3
include/openmc/source.h
View file

@ -58,6 +58,9 @@ extern "C" void initialize_source();
//! \return Sampled source site
extern "C" Bank sample_external_source();
//! Fill source bank at end of generation for fixed source simulations
void fill_source_bank_fixedsource();
} // namespace openmc
#endif // OPENMC_SOURCE_H

3
include/openmc/state_point.h
View file

@ -9,10 +9,11 @@
namespace openmc {
void write_source_point(const char* filename);
extern "C" void write_source_bank(hid_t group_id, Bank* source_bank);
extern "C" void read_source_bank(hid_t group_id, Bank* source_bank);
extern "C" void write_tally_results_nr(hid_t file_id);
extern "C" void restart_set_keff();
} // namespace openmc
#endif // OPENMC_STATE_POINT_H

7
include/openmc/tallies/tally.h
View file

@ -14,9 +14,12 @@ namespace openmc {
extern "C" double total_weight;
// Threadprivate variables
extern "C" double global_tally_absorption;
#pragma omp threadprivate(global_tally_absorption)
extern "C" double global_tally_collision;
extern "C" double global_tally_tracklength;
extern "C" double global_tally_leakage;
#pragma omp threadprivate(global_tally_absorption, global_tally_collision, \
global_tally_tracklength, global_tally_leakage)
//==============================================================================
// Non-member functions

13
include/openmc/timer.h
View file

@ -38,9 +38,22 @@ private:
// Global variables
//==============================================================================
extern Timer time_active;
extern Timer time_bank;
extern Timer time_bank_sample;
extern Timer time_bank_sendrecv;
extern Timer time_finalize;
extern Timer time_inactive;
extern Timer time_initialize;
extern Timer time_tallies;
extern Timer time_total;
extern Timer time_transport;
//==============================================================================
// Non-member functions
//==============================================================================
void reset_timers();
} // namespace openmc

240
src/api.F90
View file

@ -2,183 +2,20 @@ module openmc_api
use, intrinsic :: ISO_C_BINDING
use bank_header, only: openmc_source_bank
use constants, only: K_BOLTZMANN
use eigenvalue, only: k_sum, openmc_get_keff
use constants
use error
use geometry, only: find_cell
use geometry_header
use hdf5_interface
use material_header
use math
use message_passing
use nuclide_header
use initialize, only: openmc_init_f
use particle_header
use plot, only: openmc_plot_geometry
use random_lcg, only: openmc_get_seed, openmc_set_seed
use settings
use simulation_header
use state_point, only: openmc_statepoint_write
use tally_header
use tally_filter_header
use tally_filter
use tally, only: openmc_tally_allocate
use simulation
use string, only: to_f_string
use timer_header
use volume_calc, only: openmc_calculate_volumes
use string, only: to_str
#ifdef DAGMC
use dagmc_header, only: free_memory_dagmc
#endif
implicit none
private
public :: openmc_calculate_volumes
public :: openmc_cell_get_id
public :: openmc_cell_set_id
public :: openmc_extend_filters
public :: openmc_extend_cells
public :: openmc_extend_materials
public :: openmc_extend_tallies
public :: openmc_filter_get_id
public :: openmc_filter_get_type
public :: openmc_filter_set_id
public :: openmc_filter_set_type
public :: openmc_finalize
public :: openmc_find_cell
public :: openmc_get_cell_index
public :: openmc_get_keff
public :: openmc_get_filter_index
public :: openmc_get_filter_next_id
public :: openmc_get_material_index
public :: openmc_get_nuclide_index
public :: openmc_get_seed
public :: openmc_get_tally_index
public :: openmc_get_tally_next_id
public :: openmc_global_tallies
public :: openmc_hard_reset
public :: openmc_init_f
public :: openmc_load_nuclide
public :: openmc_material_add_nuclide
public :: openmc_material_get_id
public :: openmc_material_get_densities
public :: openmc_material_set_density
public :: openmc_material_set_densities
public :: openmc_material_set_id
public :: openmc_next_batch
public :: openmc_nuclide_name
public :: openmc_plot_geometry
public :: openmc_reset
public :: openmc_set_seed
public :: openmc_simulation_finalize
public :: openmc_simulation_init
public :: openmc_source_bank
public :: openmc_tally_allocate
public :: openmc_tally_get_estimator
public :: openmc_tally_get_id
public :: openmc_tally_get_filters
public :: openmc_tally_get_n_realizations
public :: openmc_tally_get_nuclides
public :: openmc_tally_get_scores
public :: openmc_tally_get_type
public :: openmc_tally_results
public :: openmc_tally_set_estimator
public :: openmc_tally_set_filters
public :: openmc_tally_set_id
public :: openmc_tally_set_nuclides
public :: openmc_tally_set_scores
public :: openmc_tally_set_type
contains
!===============================================================================
! OPENMC_FINALIZE frees up memory by deallocating arrays and resetting global
! variables
!===============================================================================
function openmc_finalize() result(err) bind(C)
integer(C_INT) :: err
interface
subroutine openmc_free_bank() bind(C)
end subroutine openmc_free_bank
end interface
! Clear results
err = openmc_reset()
! Reset global variables
assume_separate = .false.
check_overlaps = .false.
confidence_intervals = .false.
create_fission_neutrons = .true.
electron_treatment = ELECTRON_LED
energy_cutoff(:) = [ZERO, 1000.0_8, ZERO, ZERO]
energy_max(:) = [INFINITY, INFINITY]
energy_min(:) = [ZERO, ZERO]
entropy_on = .false.
gen_per_batch = 1
index_entropy_mesh = -1
index_ufs_mesh = -1
keff = ONE
legendre_to_tabular = .true.
legendre_to_tabular_points = C_NONE
n_batch_interval = 1
n_lost_particles = 0
n_particles = -1
n_source_points = 0
n_state_points = 0
n_tallies = 0
output_summary = .true.
output_tallies = .true.
particle_restart_run = .false.
photon_transport = .false.
pred_batches = .false.
reduce_tallies = .true.
res_scat_on = .false.
res_scat_method = RES_SCAT_ARES
res_scat_energy_min = 0.01_8
res_scat_energy_max = 1000.0_8
restart_run = .false.
root_universe = -1
run_CE = .true.
run_mode = -1
dagmc = .false.
satisfy_triggers = .false.
call openmc_set_seed(DEFAULT_SEED)
source_latest = .false.
source_separate = .false.
source_write = .true.
survival_biasing = .false.
temperature_default = 293.6_8
temperature_method = TEMPERATURE_NEAREST
temperature_multipole = .false.
temperature_range = [ZERO, ZERO]
temperature_tolerance = 10.0_8
total_gen = 0
trigger_on = .false.
ufs = .false.
urr_ptables_on = .true.
verbosity = 7
weight_cutoff = 0.25_8
weight_survive = ONE
write_all_tracks = .false.
write_initial_source = .false.
! Deallocate arrays
call free_memory()
err = 0
#ifdef OPENMC_MPI
! Free all MPI types
call openmc_free_bank()
#endif
end function openmc_finalize
!===============================================================================
! OPENMC_FIND_CELL determines what cell contains a given point in space
!===============================================================================
@ -214,85 +51,26 @@ contains
end function openmc_find_cell
!===============================================================================
! OPENMC_HARD_RESET reset tallies and timers as well as the pseudorandom
! generator state
!===============================================================================
function openmc_hard_reset() result(err) bind(C)
integer(C_INT) :: err
! Reset all tallies and timers
err = openmc_reset()
! Reset total generations and keff guess
keff = ONE
total_gen = 0
! Reset the random number generator state
call openmc_set_seed(DEFAULT_SEED)
end function openmc_hard_reset
!===============================================================================
! OPENMC_RESET resets tallies and timers
!===============================================================================
function openmc_reset() result(err) bind(C)
integer(C_INT) :: err
integer :: i
if (allocated(tallies)) then
do i = 1, size(tallies)
associate (t => tallies(i) % obj)
t % n_realizations = 0
if (allocated(t % results)) then
t % results(:, :, :) = ZERO
end if
end associate
end do
end if
! Reset global tallies
n_realizations = 0
if (allocated(global_tallies)) then
global_tallies(:, :) = ZERO
end if
k_col_abs = ZERO
k_col_tra = ZERO
k_abs_tra = ZERO
k_sum(:) = ZERO
! Reset timers
call time_total % reset()
call time_total % reset()
call time_initialize % reset()
call time_read_xs % reset()
call time_unionize % reset()
call time_tallies % reset()
call time_inactive % reset()
call time_active % reset()
call time_transport % reset()
call time_finalize % reset()
call reset_timers()
err = 0
end function openmc_reset
!===============================================================================
! FREE_MEMORY deallocates and clears all global allocatable arrays in the
! program
!===============================================================================
subroutine free_memory()
subroutine free_memory() bind(C)
use bank_header
use cmfd_header
use geometry_header
use material_header
use photon_header
use plot_header
use sab_header
use settings
use simulation_header
use surface_header
use tally_derivative_header
use tally_filter_header
use tally_header
use trigger_header
use volume_header

6
src/cmfd_execute.F90
View file

@ -29,7 +29,7 @@ contains
! EXECUTE_CMFD runs the CMFD calculation
!==============================================================================
subroutine execute_cmfd()
subroutine execute_cmfd() bind(C)
use cmfd_data, only: set_up_cmfd
use cmfd_solver, only: cmfd_solver_execute
@ -73,7 +73,7 @@ contains
! CMFD_INIT_BATCH handles cmfd options at the start of every batch
!==============================================================================
subroutine cmfd_init_batch()
subroutine cmfd_init_batch() bind(C)
! Check to activate CMFD diffusion and possible feedback
! this guarantees that when cmfd begins at least one batch of tallies are
@ -375,7 +375,7 @@ contains
! CMFD_TALLY_INIT
!===============================================================================
subroutine cmfd_tally_init()
subroutine cmfd_tally_init() bind(C)
integer :: i
if (cmfd_run) then
do i = 1, size(cmfd_tallies)

2
src/cmfd_header.F90
View file

@ -143,7 +143,7 @@ module cmfd_header
logical, public :: cmfd_run_adjoint = .false.
! CMFD run logicals
logical, public :: cmfd_on = .false.
logical(C_BOOL), public, bind(C) :: cmfd_on = .false.
! CMFD display info
character(len=25), public :: cmfd_display = 'balance'

234
src/eigenvalue.F90
View file

@ -2,241 +2,21 @@ module eigenvalue
use, intrinsic :: ISO_C_BINDING
use algorithm, only: binary_search
use constants, only: ZERO
use error, only: fatal_error, warning
use math, only: t_percentile
use message_passing
use random_lcg, only: prn, set_particle_seed, advance_prn_seed
use settings
use bank_header
use simulation_header
use string, only: to_str
use tally_header
use timer_header
implicit none
real(8) :: k_sum(2) = ZERO ! Used to reduce sum and sum_sq
interface
subroutine calculate_generation_keff() bind(C)
end subroutine
function openmc_get_keff(k_combined) result(err) bind(C)
import C_INT, C_DOUBLE
real(C_DOUBLE), intent(out) :: k_combined(2)
integer(C_INT) :: err
end function
end interface
contains
!===============================================================================
! CALCULATE_AVERAGE_KEFF calculates the mean and standard deviation of the mean
! of k-effective during active generations and broadcasts the mean to all
! processors
!===============================================================================
subroutine calculate_average_keff()
integer :: i ! overall generation within simulation
integer :: n ! number of active generations
real(8) :: alpha ! significance level for CI
real(8) :: t_value ! t-value for confidence intervals
! Determine overall generation and number of active generations
i = overall_generation()
if (current_batch > n_inactive) then
n = gen_per_batch*n_realizations + current_gen
else
n = 0
end if
if (n <= 0) then
! For inactive generations, use current generation k as estimate for next
! generation
keff = k_generation(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
! Determine mean
keff = k_sum(1) / n
if (n > 1) then
if (confidence_intervals) then
! Calculate t-value for confidence intervals
alpha = ONE - CONFIDENCE_LEVEL
t_value = t_percentile(ONE - alpha/TWO, n - 1)
else
t_value = ONE
end if
! Standard deviation of the sample mean of k
keff_std = t_value * sqrt((k_sum(2)/n - keff**2) / (n - 1))
end if
end if
end subroutine calculate_average_keff
!===============================================================================
! OPENMC_GET_KEFF calculates a minimum variance estimate of k-effective based on
! a linear combination of the collision, absorption, and tracklength
! estimates. The theory behind this can be found in M. Halperin, "Almost
! linearly-optimum combination of unbiased estimates," J. Am. Stat. Assoc., 56,
! 36-43 (1961), doi:10.1080/01621459.1961.10482088. The implementation here
! follows that described in T. Urbatsch et al., "Estimation and interpretation
! of keff confidence intervals in MCNP," Nucl. Technol., 111, 169-182 (1995).
!===============================================================================
function openmc_get_keff(k_combined) result(err) bind(C)
real(C_DOUBLE), intent(out) :: k_combined(2)
integer(C_INT) :: err
integer :: l ! loop index
integer :: i, j, k ! indices referring to collision, absorption, or track
real(8) :: n ! number of realizations
real(8) :: kv(3) ! vector of k-effective estimates
real(8) :: cov(3,3) ! sample covariance matrix
real(8) :: f ! weighting factor
real(8) :: g ! sum of weighting factors
real(8) :: S(3) ! sums used for variance calculation
k_combined = ZERO
! Make sure we have at least four realizations. Notice that at the end,
! there is a N-3 term in a denominator.
if (n_realizations <= 3) then
err = -1
return
end if
! Initialize variables
n = real(n_realizations, 8)
! Copy estimates of k-effective and its variance (not variance of the mean)
kv(1) = global_tallies(RESULT_SUM, K_COLLISION) / n
kv(2) = global_tallies(RESULT_SUM, K_ABSORPTION) / n
kv(3) = global_tallies(RESULT_SUM, K_TRACKLENGTH) / n
cov(1, 1) = (global_tallies(RESULT_SUM_SQ, K_COLLISION) - &
n * kv(1) * kv(1)) / (n - ONE)
cov(2, 2) = (global_tallies(RESULT_SUM_SQ, K_ABSORPTION) - &
n * kv(2) * kv(2)) / (n - ONE)
cov(3, 3) = (global_tallies(RESULT_SUM_SQ, K_TRACKLENGTH) - &
n * kv(3) * kv(3)) / (n - ONE)
! Calculate covariances based on sums with Bessel's correction
cov(1, 2) = (k_col_abs - n * kv(1) * kv(2)) / (n - ONE)
cov(1, 3) = (k_col_tra - n * kv(1) * kv(3)) / (n - ONE)
cov(2, 3) = (k_abs_tra - n * kv(2) * kv(3)) / (n - ONE)
cov(2, 1) = cov(1, 2)
cov(3, 1) = cov(1, 3)
cov(3, 2) = cov(2, 3)
! Check to see if two estimators are the same; this is guaranteed to happen
! in MG-mode with survival biasing when the collision and absorption
! estimators are the same, but can theoretically happen at anytime.
! If it does, the standard estimators will produce floating-point
! exceptions and an expression specifically derived for the combination of
! two estimators (vice three) should be used instead.
! First we will identify if there are any matching estimators
if ((abs(kv(1) - kv(2)) / kv(1) < FP_REL_PRECISION) .and. &
(abs(cov(1, 1) - cov(2, 2)) / cov(1, 1) < FP_REL_PRECISION)) then
! 1 and 2 match, so only use 1 and 3 in our comparisons
i = 1
j = 3
else if ((abs(kv(1) - kv(3)) / kv(1) < FP_REL_PRECISION) .and. &
(abs(cov(1, 1) - cov(3, 3)) / cov(1, 1) < FP_REL_PRECISION)) then
! 1 and 3 match, so only use 1 and 2 in our comparisons
i = 1
j = 2
else if ((abs(kv(2) - kv(3)) / kv(2) < FP_REL_PRECISION) .and. &
(abs(cov(2, 2) - cov(3, 3)) / cov(2, 2) < FP_REL_PRECISION)) then
! 2 and 3 match, so only use 1 and 2 in our comparisons
i = 1
j = 2
else
! No two estimators match, so set i to 0 and this will be the indicator
! to use all three estimators.
i = 0
end if
if (i == 0) then
! Use three estimators as derived in the paper by Urbatsch
! Initialize variables
g = ZERO
S = ZERO
do l = 1, 3
! Permutations of estimates
if (l == 1) then
! i = collision, j = absorption, k = tracklength
i = 1
j = 2
k = 3
elseif (l == 2) then
! i = absortion, j = tracklength, k = collision
i = 2
j = 3
k = 1
elseif (l == 3) then
! i = tracklength, j = collision, k = absorption
i = 3
j = 1
k = 2
end if
! Calculate weighting
f = cov(j, j) * (cov(k, k) - cov(i, k)) - cov(k, k) * cov(i, j) + &
cov(j, k) * (cov(i, j) + cov(i, k) - cov(j, k))
! Add to S sums for variance of combined estimate
S(1) = S(1) + f * cov(1, l)
S(2) = S(2) + (cov(j, j) + cov(k, k) - TWO * cov(j, k)) * kv(l) * kv(l)
S(3) = S(3) + (cov(k, k) + cov(i, j) - cov(j, k) - &
cov(i, k)) * kv(l) * kv(j)
! Add to sum for combined k-effective
k_combined(1) = k_combined(1) + f * kv(l)
g = g + f
end do
! Complete calculations of S sums
S = (n - ONE) * S
S(1) = (n - ONE)**2 * S(1)
! Calculate combined estimate of k-effective
k_combined(1) = k_combined(1) / g
! Calculate standard deviation of combined estimate
g = (n - ONE)**2 * g
k_combined(2) = sqrt(S(1) / &
(g * n * (n - THREE)) * (ONE + n * ((S(2) - TWO * S(3)) / g)))
else
! Use only two estimators
! These equations are derived analogously to that done in the paper by
! Urbatsch, but are simpler than for the three estimators case since the
! block matrices of the three estimator equations reduces to scalars here
! Store the commonly used term
f = kv(i) - kv(j)
g = cov(i, i) + cov(j, j) - TWO * cov(i, j)
! Calculate combined estimate of k-effective
k_combined(1) = kv(i) - (cov(i, i) - cov(i, j)) / g * f
! Calculate standard deviation of combined estimate
k_combined(2) = (cov(i, i) * cov(j, j) - cov(i, j) * cov(i, j)) * &
(g + n * f * f) / (n * (n - TWO) * g * g)
k_combined(2) = sqrt(k_combined(2))
end if
err = 0
end function openmc_get_keff
#ifdef _OPENMP
!===============================================================================
! JOIN_BANK_FROM_THREADS joins threadprivate fission banks into a single fission
@ -244,7 +24,7 @@ contains
! to preserve the order of the bank when using varying numbers of threads.
!===============================================================================
subroutine join_bank_from_threads()
subroutine join_bank_from_threads() bind(C)
integer(8) :: total ! total number of fission bank sites
integer :: i ! loop index for threads

201
src/eigenvalue.cpp
View file

@ -1,5 +1,6 @@
#include "openmc/eigenvalue.h"
#include "xtensor/xbuilder.hpp"
#include "xtensor/xmath.hpp"
#include "xtensor/xtensor.hpp"
#include "xtensor/xview.hpp"
@ -8,6 +9,7 @@
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/hdf5_interface.h"
#include "openmc/math_functions.h"
#include "openmc/mesh.h"
#include "openmc/message_passing.h"
#include "openmc/random_lcg.h"
@ -18,6 +20,8 @@
#include "openmc/tallies/tally.h"
#include <algorithm> // for min
#include <array>
#include <cmath> // for sqrt, abs, pow
#include <string>
namespace openmc {
@ -27,6 +31,7 @@ namespace openmc {
//==============================================================================
double keff_generation;
std::array<double, 2> k_sum;
std::vector<double> entropy;
xt::xtensor<double, 1> source_frac;
@ -129,7 +134,7 @@ void synchronize_bank()
// Allocate temporary source bank
int64_t index_temp = 0;
Bank temp_sites[3*simulation::work];
std::vector<Bank> temp_sites(3*simulation::work);
for (int64_t i = 0; i < n_bank; ++i) {
// If there are less than n_particles particles banked, automatically add
@ -291,13 +296,200 @@ void synchronize_bank()
MPI_Waitall(n_request, requests.data(), MPI_STATUSES_IGNORE);
#else
std::copy(temp_sites, temp_sites + settings::n_particles, source_bank);
std::copy(temp_sites.data(), temp_sites.data() + settings::n_particles, source_bank);
#endif
time_bank_sendrecv.stop();
time_bank.stop();
}
void calculate_average_keff()
{
// Determine overall generation and number of active generations
int i = overall_generation() - 1;
int n;
if (simulation::current_batch > settings::n_inactive) {
n = settings::gen_per_batch*n_realizations + simulation::current_gen;
} else {
n = 0;
}
if (n <= 0) {
// For inactive generations, use current generation k as estimate for next
// generation
simulation::keff = simulation::k_generation[i];
} else {
// Sample mean of keff
k_sum[0] += simulation::k_generation[i];
k_sum[1] += std::pow(simulation::k_generation[i], 2);
// Determine mean
simulation::keff = k_sum[0] / n;
if (n > 1) {
double t_value;
if (settings::confidence_intervals) {
// Calculate t-value for confidence intervals
double alpha = 1.0 - CONFIDENCE_LEVEL;
t_value = t_percentile_c(1.0 - alpha/2.0, n - 1);
} else {
t_value = 1.0;
}
// Standard deviation of the sample mean of k
simulation::keff_std = t_value * std::sqrt((k_sum[1]/n -
std::pow(simulation::keff, 2)) / (n - 1));
}
}
}
int openmc_get_keff(double* k_combined)
{
k_combined[0] = 0.0;
k_combined[1] = 0.0;
// Make sure we have at least four realizations. Notice that at the end,
// there is a N-3 term in a denominator.
if (n_realizations <= 3) {
return -1;
}
// Initialize variables
int64_t n = n_realizations;
// Copy estimates of k-effective and its variance (not variance of the mean)
auto gt = global_tallies();
std::array<double, 3> kv {};
xt::xtensor<double, 2> cov = xt::zeros<double>({3, 3});
kv[0] = gt(K_COLLISION, RESULT_SUM) / n;
kv[1] = gt(K_ABSORPTION, RESULT_SUM) / n;
kv[2] = gt(K_TRACKLENGTH, RESULT_SUM) / n;
cov(0, 0) = (gt(K_COLLISION, RESULT_SUM_SQ) - n*kv[0]*kv[0]) / (n - 1);
cov(1, 1) = (gt(K_ABSORPTION, RESULT_SUM_SQ) - n*kv[1]*kv[1]) / (n - 1);
cov(2, 2) = (gt(K_TRACKLENGTH, RESULT_SUM_SQ) - n*kv[2]*kv[2]) / (n - 1);
// Calculate covariances based on sums with Bessel's correction
cov(0, 1) = (simulation::k_col_abs - n * kv[0] * kv[1]) / (n - 1);
cov(0, 2) = (simulation::k_col_tra - n * kv[0] * kv[2]) / (n - 1);
cov(1, 2) = (simulation::k_abs_tra - n * kv[1] * kv[2]) / (n - 1);
cov(1, 0) = cov(0, 1);
cov(2, 0) = cov(0, 2);
cov(2, 1) = cov(1, 2);
// Check to see if two estimators are the same; this is guaranteed to happen
// in MG-mode with survival biasing when the collision and absorption
// estimators are the same, but can theoretically happen at anytime.
// If it does, the standard estimators will produce floating-point
// exceptions and an expression specifically derived for the combination of
// two estimators (vice three) should be used instead.
// First we will identify if there are any matching estimators
int i, j, k;
if ((std::abs(kv[0] - kv[1]) / kv[0] < FP_REL_PRECISION) &&
(std::abs(cov(0, 0) - cov(1, 1)) / cov(0, 0) < FP_REL_PRECISION)) {
// 0 and 1 match, so only use 0 and 2 in our comparisons
i = 0;
j = 2;
} else if ((std::abs(kv[0] - kv[2]) / kv[0] < FP_REL_PRECISION) &&
(std::abs(cov(0, 0) - cov(2, 2)) / cov(0, 0) < FP_REL_PRECISION)) {
// 0 and 2 match, so only use 0 and 1 in our comparisons
i = 0;
j = 1;
} else if ((std::abs(kv[1] - kv[2]) / kv[1] < FP_REL_PRECISION) &&
(std::abs(cov(1, 1) - cov(2, 2)) / cov(1, 1) < FP_REL_PRECISION)) {
// 1 and 2 match, so only use 0 and 1 in our comparisons
i = 0;
j = 1;
} else {
// No two estimators match, so set i to -1 and this will be the indicator
// to use all three estimators.
i = -1;
}
if (i == -1) {
// Use three estimators as derived in the paper by Urbatsch
// Initialize variables
double g = 0.0;
std::array<double, 3> S {};
for (int l = 0; l < 3; ++l) {
// Permutations of estimates
switch (l) {
case 0:
// i = collision, j = absorption, k = tracklength
i = 0;
j = 1;
k = 2;
break;
case 1:
// i = absortion, j = tracklength, k = collision
i = 1;
j = 2;
k = 0;
break;
case 2:
// i = tracklength, j = collision, k = absorption
i = 2;
j = 0;
k = 1;
break;
}
// Calculate weighting
double f = cov(j, j) * (cov(k, k) - cov(i, k)) - cov(k, k) * cov(i, j) +
cov(j, k) * (cov(i, j) + cov(i, k) - cov(j, k));
// Add to S sums for variance of combined estimate
S[0] += f * cov(0, l);
S[1] += (cov(j, j) + cov(k, k) - 2.0 * cov(j, k)) * kv[l] * kv[l];
S[2] += (cov(k, k) + cov(i, j) - cov(j, k) - cov(i, k)) * kv[l] * kv[j];
// Add to sum for combined k-effective
k_combined[0] += f * kv[l];
g += f;
}
// Complete calculations of S sums
for (auto& S_i : S) {
S_i *= (n - 1);
}
S[0] *= (n - 1)*(n - 1);
// Calculate combined estimate of k-effective
k_combined[0] /= g;
// Calculate standard deviation of combined estimate
g *= (n - 1)*(n - 1);
k_combined[1] = std::sqrt(S[0] / (g*n*(n - 3)) *
(1 + n*((S[1] - 2*S[2]) / g)));
} else {
// Use only two estimators
// These equations are derived analogously to that done in the paper by
// Urbatsch, but are simpler than for the three estimators case since the
// block matrices of the three estimator equations reduces to scalars here
// Store the commonly used term
double f = kv[i] - kv[j];
double g = cov(i, i) + cov(j, j) - 2.0*cov(i, j);
// Calculate combined estimate of k-effective
k_combined[0] = kv[i] - (cov(i, i) - cov(i, j)) / g * f;
// Calculate standard deviation of combined estimate
k_combined[1] = (cov(i, i)*cov(j, j) - cov(i, j)*cov(i, j)) *
(g + n*f*f) / (n*(n - 2)*g*g);
k_combined[1] = std::sqrt(k_combined[1]);
}
return 0;
}
void shannon_entropy()
{
// Get pointer to entropy mesh
@ -430,6 +622,10 @@ extern "C" void read_eigenvalue_hdf5(hid_t group)
read_dataset(group, "k_abs_tra", simulation::k_abs_tra);
}
//==============================================================================
// Fortran compatibility
//==============================================================================
extern "C" double entropy_c(int i)
{
return entropy.at(i - 1);
@ -440,5 +636,6 @@ extern "C" void entropy_clear()
entropy.clear();
}
extern "C" void k_sum_reset() { k_sum.fill(0.0); }
} // namespace openmc

129
src/finalize.cpp
View file

@ -1,14 +1,137 @@
#include "openmc/finalize.h"
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/eigenvalue.h"
#include "openmc/geometry.h"
#include "openmc/message_passing.h"
#include "openmc/nuclide.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/timer.h"
#include "openmc/tallies/tally.h"
namespace openmc {
using namespace openmc;
void openmc_free_bank()
// Functions defined in Fortran
extern "C" void free_memory();
extern "C" void reset_timers_f();
int openmc_finalize()
{
// Clear results
openmc_reset();
// Reset global variables
settings::assume_separate = false;
settings::check_overlaps = false;
settings::confidence_intervals = false;
settings::create_fission_neutrons = true;
settings::electron_treatment = ELECTRON_LED;
settings::energy_cutoff = {0.0, 1000.0, 0.0, 0.0};
settings::entropy_on = false;
settings::gen_per_batch = 1;
settings::index_entropy_mesh = -1;
settings::index_ufs_mesh = -1;
settings::legendre_to_tabular = true;
settings::legendre_to_tabular_points = -1;
settings::n_particles = -1;
settings::output_summary = true;
settings::output_tallies = true;
settings::particle_restart_run = false;
settings::photon_transport = false;
settings::reduce_tallies = true;
settings::res_scat_on = false;
settings::res_scat_method = RES_SCAT_ARES;
settings::res_scat_energy_min = 0.01;
settings::res_scat_energy_max = 1000.0;
settings::restart_run = false;
settings::run_CE = true;
settings::run_mode = -1;
settings::dagmc = false;
settings::source_latest = false;
settings::source_separate = false;
settings::source_write = true;
settings::survival_biasing = false;
settings::temperature_default = 293.6;
settings::temperature_method = TEMPERATURE_NEAREST;
settings::temperature_multipole = false;
settings::temperature_range = {0.0, 0.0};
settings::temperature_tolerance = 10.0;
settings::trigger_on = false;
settings::trigger_predict = false;
settings::trigger_batch_interval = 1;
settings::ufs_on = false;
settings::urr_ptables_on = true;
settings::verbosity = 7;
settings::weight_cutoff = 0.25;
settings::weight_survive = 1.0;
settings::write_all_tracks = false;
settings::write_initial_source = false;
simulation::keff = 1.0;
simulation::n_lost_particles = 0;
simulation::satisfy_triggers = false;
simulation::total_gen = 0;
energy_max = {INFTY, INFTY};
energy_min = {0.0, 0.0};
n_tallies = 0;
openmc_root_universe = -1;
openmc_set_seed(DEFAULT_SEED);
// Deallocate arrays
free_memory();
// Free all MPI types
#ifdef OPENMC_MPI
MPI_Type_free(&mpi::bank);
#endif
return 0;
}
} // namespace openmc
int openmc_reset()
{
for (int i = 1; i <= n_tallies; ++i) {
openmc_tally_reset(i);
}
// Reset global tallies (can't really use global_tallies() right now because
// it doesn't have any information about whether the underlying buffer was
// allocated)
n_realizations = 0;
double* buffer = nullptr;
openmc_global_tallies(&buffer);
if (buffer) {
for (int i = 0; i < 3*N_GLOBAL_TALLIES; ++i) {
buffer[i] = 0.0;
}
}
simulation::k_col_abs = 0.0;
simulation::k_col_tra = 0.0;
simulation::k_abs_tra = 0.0;
k_sum = {0.0, 0.0};
// Reset timers
reset_timers();
reset_timers_f();
return 0;
}
int openmc_hard_reset()
{
// Reset all tallies and timers
openmc_reset();
// Reset total generations and keff guess
simulation::keff = 1.0;
simulation::total_gen = 0;
// Reset the random number generator state
openmc_set_seed(DEFAULT_SEED);
return 0;
}

21
src/hdf5_interface.cpp
View file

@ -148,13 +148,22 @@ dataset_typesize(hid_t dset)
void
ensure_exists(hid_t group_id, const char* name)
ensure_exists(hid_t obj_id, const char* name, bool attribute)
{
if (!object_exists(group_id, name)) {
std::stringstream err_msg;
err_msg << "Object \"" << name << "\" does not exist in group "
<< object_name(group_id);
fatal_error(err_msg);
if (attribute) {
if (!attribute_exists(obj_id, name)) {
std::stringstream err_msg;
err_msg << "Attribute \"" << name << "\" does not exist in object "
<< object_name(obj_id);
fatal_error(err_msg);
}
} else {
if (!object_exists(obj_id, name)) {
std::stringstream err_msg;
err_msg << "Object \"" << name << "\" does not exist in object "
<< object_name(obj_id);
fatal_error(err_msg);
}
}
}

59
src/initialize.F90
View file

@ -2,18 +2,8 @@ module initialize
use, intrinsic :: ISO_C_BINDING
#ifdef _OPENMP
use omp_lib
#endif
use bank_header, only: Bank
use constants
use input_xml, only: read_input_xml
use message_passing
use random_lcg, only: openmc_set_seed
use settings
use string, only: ends_with, to_f_string
use timer_header
use string, only: to_f_string
implicit none
@ -42,56 +32,11 @@ module initialize
contains
!===============================================================================
! OPENMC_INIT takes care of all initialization tasks, i.e. reading
! from command line, reading xml input files, initializing random
! number seeds, reading cross sections, initializing starting source,
! setting up timers, etc.
!===============================================================================
function openmc_init_f() result(err) bind(C)
integer(C_INT) :: err
#ifdef _OPENMP
character(MAX_WORD_LEN) :: envvar
#endif
! Start total and initialization timer
call time_total%start()
call time_initialize%start()
#ifdef _OPENMP
! Change schedule of main parallel-do loop if OMP_SCHEDULE is set
call get_environment_variable("OMP_SCHEDULE", envvar)
if (len_trim(envvar) == 0) then
call omp_set_schedule(omp_sched_static, 0)
end if
#endif
! Read command line arguments
call read_command_line()
! Initialize random number generator -- if the user specifies a seed, it
! will be re-initialized later
call openmc_set_seed(DEFAULT_SEED)
! Read XML input files
call read_input_xml()
! Check for particle restart run
if (particle_restart_run) run_mode = MODE_PARTICLE
! Stop initialization timer
call time_initialize%stop()
err = 0
end function openmc_init_f
!===============================================================================
! READ_COMMAND_LINE reads all parameters from the command line
!===============================================================================
subroutine read_command_line()
subroutine read_command_line() bind(C)
! Arguments were already read on C++ side (initialize.cpp). Here we just
! convert the C-style strings to Fortran style

43
src/initialize.cpp
View file

@ -1,12 +1,13 @@
#include "openmc/initialize.h"
#include <cstddef>
#include <cstdlib> // for getenv
#include <cstring>
#include <sstream>
#include <string>
#ifdef _OPENMP
#include "omp.h"
#include <omp.h>
#endif
#include "openmc/capi.h"
@ -14,14 +15,17 @@
#include "openmc/error.h"
#include "openmc/hdf5_interface.h"
#include "openmc/message_passing.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/string_utils.h"
#include "openmc/timer.h"
// data/functions from Fortran side
extern "C" void openmc_init_f();
extern "C" void print_usage();
extern "C" void print_version();
extern "C" void read_command_line();
extern "C" void read_input_xml();
// Paths to various files
extern "C" {
@ -30,6 +34,8 @@ extern "C" {
int openmc_init(int argc, char* argv[], const void* intracomm)
{
using namespace openmc;
#ifdef OPENMC_MPI
// Check if intracomm was passed
MPI_Comm comm;
@ -40,15 +46,40 @@ int openmc_init(int argc, char* argv[], const void* intracomm)
}
// Initialize MPI for C++
openmc::initialize_mpi(comm);
initialize_mpi(comm);
#endif
// Parse command-line arguments
int err = openmc::parse_command_line(argc, argv);
int err = parse_command_line(argc, argv);
if (err) return err;
// Continue with rest of initialization
openmc_init_f();
// Start total and initialization timer
time_total.start();
time_initialize.start();
#ifdef _OPENMP
// If OMP_SCHEDULE is not set, default to a static schedule
char* envvar = std::getenv("OMP_SCHEDULE");
if (!envvar) {
omp_set_schedule(omp_sched_static, 0);
}
#endif
// Read command line arguments
read_command_line();
// Initialize random number generator -- if the user specifies a seed, it
// will be re-initialized later
openmc_set_seed(DEFAULT_SEED);
// Read XML input files
read_input_xml();
// Check for particle restart run
if (settings::particle_restart_run) settings::run_mode = RUN_MODE_PARTICLE;
// Stop initialization timer
time_initialize.stop();
return 0;
}

128
src/input_xml.F90
View file

@ -28,7 +28,7 @@ module input_xml
use plot_header
use random_lcg, only: prn
use surface_header
use set_header, only: SetChar
use set_header, only: SetChar, SetInt
use settings
use stl_vector, only: VectorInt, VectorReal, VectorChar
use string, only: to_lower, to_str, str_to_int, str_to_real, &
@ -112,7 +112,7 @@ contains
! geometry, materials, and tallies.
!===============================================================================
subroutine read_input_xml()
subroutine read_input_xml() bind(C)
type(VectorReal), allocatable :: nuc_temps(:) ! List of T to read for each nuclide
type(VectorReal), allocatable :: sab_temps(:) ! List of T to read for each S(a,b)
@ -207,10 +207,7 @@ contains
integer :: i
integer :: n
integer, allocatable :: temp_int_array(:)
integer :: n_tracks
type(XMLNode) :: root
type(XMLNode) :: node_sp
type(XMLNode) :: node_res_scat
type(XMLNode) :: node_vol
type(XMLNode), allocatable :: node_vol_list(:)
@ -218,113 +215,6 @@ contains
! Get proper XMLNode type given pointer
root % ptr = root_ptr
if (run_mode == MODE_EIGENVALUE) then
! Preallocate space for keff and entropy by generation
call k_generation_reserve(n_max_batches*gen_per_batch)
end if
! Particle tracks
if (check_for_node(root, "track")) then
! Make sure that there are three values per particle
n_tracks = node_word_count(root, "track")
if (mod(n_tracks, 3) /= 0) then
call fatal_error("Number of integers specified in 'track' is not &
&divisible by 3. Please provide 3 integers per particle to be &
&tracked.")
end if
! Allocate space and get list of tracks
allocate(temp_int_array(n_tracks))
call get_node_array(root, "track", temp_int_array)
! Reshape into track_identifiers
allocate(track_identifiers(3, n_tracks/3))
track_identifiers = reshape(temp_int_array, [3, n_tracks/3])
end if
! Check if the user has specified to write state points
if (check_for_node(root, "state_point")) then
! Get pointer to state_point node
node_sp = root % child("state_point")
! Determine number of batches at which to store state points
if (check_for_node(node_sp, "batches")) then
n_state_points = node_word_count(node_sp, "batches")
else
n_state_points = 0
end if
if (n_state_points > 0) then
! User gave specific batches to write state points
allocate(temp_int_array(n_state_points))
call get_node_array(node_sp, "batches", temp_int_array)
do i = 1, n_state_points
call statepoint_batch % add(temp_int_array(i))
end do
deallocate(temp_int_array)
else
! If neither were specified, write state point at last batch
n_state_points = 1
call statepoint_batch % add(n_batches)
end if
else
! If no <state_point> tag was present, by default write state point at
! last batch only
n_state_points = 1
call statepoint_batch % add(n_batches)
end if
! Check if the user has specified to write source points
if (check_for_node(root, "source_point")) then
! Get pointer to source_point node
node_sp = root % child("source_point")
! Determine number of batches at which to store source points
if (check_for_node(node_sp, "batches")) then
n_source_points = node_word_count(node_sp, "batches")
else
n_source_points = 0
end if
if (n_source_points > 0) then
! User gave specific batches to write source points
allocate(temp_int_array(n_source_points))
call get_node_array(node_sp, "batches", temp_int_array)
do i = 1, n_source_points
call sourcepoint_batch % add(temp_int_array(i))
end do
deallocate(temp_int_array)
else
! If neither were specified, write source points with state points
n_source_points = n_state_points
do i = 1, n_state_points
call sourcepoint_batch % add(statepoint_batch % get_item(i))
end do
end if
else
! If no <source_point> tag was present, by default we keep source bank in
! statepoint file and write it out at statepoints intervals
n_source_points = n_state_points
do i = 1, n_state_points
call sourcepoint_batch % add(statepoint_batch % get_item(i))
end do
end if
! If source is not seperate and is to be written out in the statepoint file,
! make sure that the sourcepoint batch numbers are contained in the
! statepoint list
if (.not. source_separate) then
do i = 1, n_source_points
if (.not. statepoint_batch % contains(sourcepoint_batch % &
get_item(i))) then
call fatal_error('Sourcepoint batches are not a subset&
& of statepoint batches.')
end if
end do
end if
! Resonance scattering parameters
if (check_for_node(root, "resonance_scattering")) then
node_res_scat = root % child("resonance_scattering")
@ -1322,7 +1212,11 @@ contains
! Read derivative attributes.
do i = 1, size(node_deriv_list)
!$omp parallel
!$omp critical (ReadTallyDeriv)
call tally_derivs(i) % from_xml(node_deriv_list(i))
!$omp end critical (ReadTallyDeriv)
!$omp end parallel
! Update tally derivative dictionary
call tally_deriv_dict % set(tally_derivs(i) % id, i)
@ -2714,6 +2608,13 @@ contains
integer(HID_T) :: file_id
character(len=MAX_WORD_LEN), allocatable :: names(:)
interface
subroutine read_mg_cross_sections_header_c(file_id) bind(C)
import HID_T
integer(HID_T), value :: file_id
end subroutine
end interface
! Check if MGXS Library exists
inquire(FILE=path_cross_sections, EXIST=file_exists)
if (.not. file_exists) then
@ -2761,6 +2662,9 @@ contains
energy_bin_avg(i) = HALF * (energy_bins(i) + energy_bins(i + 1))
end do
! Set up energy bins on C++ side
call read_mg_cross_sections_header_c(file_id)
! Get the minimum and maximum energies
energy_min(NEUTRON) = energy_bins(num_energy_groups + 1)
energy_max(NEUTRON) = energy_bins(1)

1
src/mgxs.cpp
View file

@ -26,7 +26,6 @@ namespace openmc {
std::vector<Mgxs> nuclides_MG;
std::vector<Mgxs> macro_xs;
//==============================================================================
// Mgxs base-class methods
//==============================================================================

28
src/mgxs_interface.cpp
View file

@ -8,6 +8,14 @@
namespace openmc {
//==============================================================================
// Global variable definitions
//==============================================================================
std::vector<double> energy_bins;
std::vector<double> energy_bin_avg;
std::vector<double> rev_energy_bins;
//==============================================================================
// Mgxs data loading interface methods
//==============================================================================
@ -83,6 +91,26 @@ create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[],
}
}
//==============================================================================
void read_mg_cross_sections_header_c(hid_t file_id)
{
ensure_exists(file_id, "energy_groups", true);
read_attribute(file_id, "energy_groups", num_energy_groups);
ensure_exists(file_id, "group structure", true);
read_attribute(file_id, "group structure", energy_bins);
// Create reverse energy bins
std::copy(energy_bins.crbegin(), energy_bins.crend(),
std::back_inserter(rev_energy_bins));
// Create average energies
for (int i = 0; i < energy_bins.size() - 1; ++i) {
energy_bin_avg.push_back(0.5*(energy_bins[i] + energy_bins[i+1]));
}
}
//==============================================================================
// Mgxs tracking/transport/tallying interface methods
//==============================================================================

40
src/output.F90
View file

@ -285,7 +285,7 @@ contains
! below them
!===============================================================================
subroutine print_columns()
subroutine print_columns() bind(C)
write(UNIT=ou, FMT='(2X,A9,3X)', ADVANCE='NO') "Bat./Gen."
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') " k "
@ -323,7 +323,7 @@ contains
! PRINT_GENERATION displays information for a generation of neutrons.
!===============================================================================
subroutine print_generation()
subroutine print_generation() bind(C)
integer :: i ! overall generation
integer :: n ! number of active generations
@ -360,7 +360,7 @@ contains
! multiplication factor as well as the average value if we're in active batches
!===============================================================================
subroutine print_batch_keff()
subroutine print_batch_keff() bind(C)
integer :: i ! overall generation
integer :: n ! number of active generations
@ -490,7 +490,7 @@ contains
! initialization, for computation, and for intergeneration synchronization.
!===============================================================================
subroutine print_runtime()
subroutine print_runtime() bind(C)
integer :: n_active
real(8) :: speed_inactive ! # of neutrons/second in inactive batches
@ -501,49 +501,49 @@ contains
call header("Timing Statistics", 6)
! display time elapsed for various sections
write(ou,100) "Total time for initialization", time_initialize % elapsed
write(ou,100) "Total time for initialization", time_initialize_elapsed()
write(ou,100) " Reading cross sections", time_read_xs % elapsed
write(ou,100) "Total time in simulation", time_inactive % elapsed + &
time_active % elapsed
write(ou,100) " Time in transport only", time_transport % elapsed
write(ou,100) "Total time in simulation", time_inactive_elapsed() + &
time_active_elapsed()
write(ou,100) " Time in transport only", time_transport_elapsed()
if (run_mode == MODE_EIGENVALUE) then
write(ou,100) " Time in inactive batches", time_inactive % elapsed
write(ou,100) " Time in inactive batches", time_inactive_elapsed()
end if
write(ou,100) " Time in active batches", time_active % elapsed
write(ou,100) " Time in active batches", time_active_elapsed()
if (run_mode == MODE_EIGENVALUE) then
write(ou,100) " Time synchronizing fission bank", time_bank_elapsed()
write(ou,100) " Sampling source sites", time_bank_sample_elapsed()
write(ou,100) " SEND/RECV source sites", time_bank_sendrecv_elapsed()
end if
write(ou,100) " Time accumulating tallies", time_tallies % elapsed
write(ou,100) " Time accumulating tallies", time_tallies_elapsed()
if (cmfd_run) write(ou,100) " Time in CMFD", time_cmfd % elapsed
if (cmfd_run) write(ou,100) " Building matrices", &
time_cmfdbuild % elapsed
if (cmfd_run) write(ou,100) " Solving matrices", &
time_cmfdsolve % elapsed
write(ou,100) "Total time for finalization", time_finalize % elapsed
write(ou,100) "Total time elapsed", time_total % elapsed
write(ou,100) "Total time for finalization", time_finalize_elapsed()
write(ou,100) "Total time elapsed", time_total_elapsed()
! Calculate particle rate in active/inactive batches
n_active = current_batch - n_inactive
if (restart_run) then
if (restart_batch < n_inactive) then
speed_inactive = real(n_particles * (n_inactive - restart_batch) * &
gen_per_batch) / time_inactive % elapsed
gen_per_batch) / time_inactive_elapsed()
speed_active = real(n_particles * n_active * gen_per_batch) / &
time_active % elapsed
time_active_elapsed()
else
speed_inactive = ZERO
speed_active = real(n_particles * (n_batches - restart_batch) * &
gen_per_batch) / time_active % elapsed
gen_per_batch) / time_active_elapsed()
end if
else
if (n_inactive > 0) then
speed_inactive = real(n_particles * n_inactive * gen_per_batch) / &
time_inactive % elapsed
time_inactive_elapsed()
end if
speed_active = real(n_particles * n_active * gen_per_batch) / &
time_active % elapsed
time_active_elapsed()
end if
! display calculation rate
@ -566,7 +566,7 @@ contains
! leakage rate.
!===============================================================================
subroutine print_results()
subroutine print_results() bind(C)
integer :: n ! number of realizations
real(8) :: alpha ! significance level for CI
@ -633,7 +633,7 @@ contains
! tallies and their standard deviations
!===============================================================================
subroutine write_tallies()
subroutine write_tallies() bind(C)
integer :: i ! index in tallies array
integer :: j ! level in tally hierarchy

3
src/output.cpp
View file

@ -41,8 +41,7 @@ header(const char* msg, int level) {
//==============================================================================
extern "C" void
print_overlap_check() {
void print_overlap_check() {
#ifdef OPENMC_MPI
std::vector<int64_t> temp(overlap_check_count);
int err = MPI_Reduce(temp.data(), overlap_check_count.data(),

10
src/particle.cpp
View file

@ -7,6 +7,7 @@
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/hdf5_interface.h"
#include "openmc/mgxs_interface.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
@ -92,7 +93,7 @@ Particle::initialize()
}
void
Particle::from_source(const Bank* src, bool run_CE, const double* energy_bin_avg)
Particle::from_source(const Bank* src)
{
// set defaults
initialize();
@ -106,7 +107,7 @@ Particle::from_source(const Bank* src, bool run_CE, const double* energy_bin_avg
std::copy(src->xyz, src->xyz + 3, last_xyz_current);
std::copy(src->xyz, src->xyz + 3, last_xyz);
std::copy(src->uvw, src->uvw + 3, last_uvw);
if (run_CE) {
if (settings::run_CE) {
E = src->E;
g = 0;
} else {
@ -212,10 +213,9 @@ void particle_create_secondary(Particle* p, const double* uvw, double E,
p->create_secondary(uvw, E, type, run_CE);
}
void particle_initialize(Particle* p) { p->initialize(); }
void particle_from_source(Particle* p, const Bank* src, bool run_CE,
const double* energy_bin_avg)
void particle_from_source(Particle* p, const Bank* src)
{
p->from_source(src, run_CE, energy_bin_avg);
p->from_source(src);
}
void particle_mark_as_lost(Particle* p, const char* message)
{

6
src/particle_header.F90
View file

@ -129,12 +129,10 @@ module particle_header
type(Particle), intent(inout) :: p
end subroutine particle_initialize
subroutine particle_from_source(p, src, run_CE, energy_bin_avg) bind(C)
import Particle, Bank, C_BOOL, C_DOUBLE
subroutine particle_from_source(p, src) bind(C)
import Particle, Bank, C_DOUBLE
type(Particle), intent(inout) :: p
type(Bank), intent(in) :: src
logical(C_BOOL), value :: run_CE
real(C_DOUBLE), intent(in) :: energy_bin_avg(*)
end subroutine particle_from_source
subroutine particle_mark_as_lost_c(p, message) bind(C, name='particle_mark_as_lost')

31
src/settings.F90
View file

@ -3,7 +3,6 @@ module settings
use, intrinsic :: ISO_C_BINDING
use constants
use set_header, only: SetInt
implicit none
@ -67,9 +66,6 @@ module settings
integer(C_INT32_T), bind(C) :: index_ufs_mesh
! Write source at end of simulation
logical(C_BOOL), bind(C) :: source_separate
logical(C_BOOL), bind(C) :: source_write
logical(C_BOOL), bind(C) :: source_latest
! Variance reduction settins
logical(C_BOOL), bind(C) :: survival_biasing
@ -99,25 +95,10 @@ module settings
! Particle tracks
logical(C_BOOL), bind(C) :: write_all_tracks
integer, allocatable :: track_identifiers(:,:)
! Particle restart run
logical(C_BOOL), bind(C) :: particle_restart_run
! Write out initial source
logical(C_BOOL), bind(C) :: write_initial_source
! Whether create fission neutrons or not. Only applied for MODE_FIXEDSOURCE
logical(C_BOOL), bind(C) :: create_fission_neutrons
! Information about state points to be written
integer :: n_state_points = 0
type(SetInt) :: statepoint_batch
! Information about source points to be written
integer :: n_source_points = 0
type(SetInt) :: sourcepoint_batch
character(MAX_FILE_LEN) :: path_input ! Path to input file
character(MAX_FILE_LEN) :: path_cross_sections = '' ! Path to cross_sections.xml
character(MAX_FILE_LEN) :: path_multipole ! Path to wmp library
@ -128,7 +109,6 @@ module settings
! Various output options
logical(C_BOOL), bind(C) :: output_summary
logical(C_BOOL), bind(C) :: output_tallies
! Resonance scattering settings
logical(C_BOOL), bind(C) :: res_scat_on ! is resonance scattering treated?
@ -150,11 +130,14 @@ contains
!===============================================================================
subroutine free_memory_settings()
if (allocated(res_scat_nuclides)) deallocate(res_scat_nuclides)
if (allocated(track_identifiers)) deallocate(track_identifiers)
interface
subroutine free_memory_settings_c() bind(C)
end subroutine
end interface
call statepoint_batch % clear()
call sourcepoint_batch % clear()
if (allocated(res_scat_nuclides)) deallocate(res_scat_nuclides)
call free_memory_settings_c()
end subroutine free_memory_settings
end module settings

81
src/settings.cpp
View file

@ -9,9 +9,11 @@
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/container_util.h"
#include "openmc/distribution.h"
#include "openmc/distribution_multi.h"
#include "openmc/distribution_spatial.h"
#include "openmc/eigenvalue.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/mesh.h"
@ -77,7 +79,7 @@ int32_t gen_per_batch {1};
int64_t n_particles {-1};
int electron_treatment {ELECTRON_TTB};
double energy_cutoff[4] {0.0, 1000.0, 0.0, 0.0};
std::array<double, 4> energy_cutoff {0.0, 1000.0, 0.0, 0.0};
int legendre_to_tabular_points {C_NONE};
int max_order {0};
int n_log_bins {8000};
@ -86,13 +88,16 @@ int res_scat_method {RES_SCAT_ARES};
double res_scat_energy_min {0.01};
double res_scat_energy_max {1000.0};
int run_mode {-1};
std::unordered_set<int> sourcepoint_batch;
std::unordered_set<int> statepoint_batch;
int temperature_method {TEMPERATURE_NEAREST};
double temperature_tolerance {10.0};
double temperature_default {293.6};
double temperature_range[2] {0.0, 0.0};
std::array<double, 2> temperature_range {0.0, 0.0};
int trace_batch;
int trace_gen;
int64_t trace_particle;
std::vector<std::array<int, 3>> track_identifiers;
int trigger_batch_interval {1};
int verbosity {7};
double weight_cutoff {0.25};
@ -141,7 +146,10 @@ void get_run_parameters(pugi::xml_node node_base)
gen_per_batch = std::stoi(get_node_value(node_base, "generations_per_batch"));
}
// TODO: Preallocate space for keff and entropy by generation
// Preallocate space for keff and entropy by generation
int m = settings::n_max_batches * settings::gen_per_batch;
simulation::k_generation.reserve(m);
entropy.reserve(m);
// Get the trigger information for keff
if (check_for_node(node_base, "keff_trigger")) {
@ -491,7 +499,7 @@ void read_settings_xml()
// Particle tracks
if (check_for_node(root, "track")) {
// Get values and make sure there are three per particle
auto temp = get_node_array<int64_t>(root, "track");
auto temp = get_node_array<int>(root, "track");
if (temp.size() % 3 != 0) {
fatal_error("Number of integers specified in 'track' is not "
"divisible by 3. Please provide 3 integers per particle to be "
@ -499,8 +507,11 @@ void read_settings_xml()
}
// Reshape into track_identifiers
//allocate(track_identifiers(3, n_tracks/3))
//track_identifiers = reshape(temp_int_array, [3, n_tracks/3])
int n_tracks = temp.size() / 3;
for (int i = 0; i < n_tracks; ++i) {
track_identifiers.push_back({temp[3*i], temp[3*i + 1],
temp[3*i + 2]});
}
}
// Read meshes
@ -539,7 +550,7 @@ void read_settings_xml()
// 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
int n = std::ceil(std::pow(settings::n_particles / 20.0, 1.0/3.0));
int n = std::ceil(std::pow(n_particles / 20.0, 1.0/3.0));
m.shape_ = {n, n, n};
m.n_dimension_ = 3;
@ -548,7 +559,7 @@ void read_settings_xml()
}
// Turn on Shannon entropy calculation
settings::entropy_on = true;
entropy_on = true;
}
// Uniform fission source weighting mesh
@ -581,18 +592,48 @@ void read_settings_xml()
if (index_ufs_mesh >= 0) {
// Turn on uniform fission source weighting
settings::ufs_on = true;
ufs_on = true;
}
// TODO: Read <state_point>
// Check if the user has specified to write state points
if (check_for_node(root, "state_point")) {
// Get pointer to state_point node
auto node_sp = root.child("state_point");
// Determine number of batches at which to store state points
if (check_for_node(node_sp, "batches")) {
// User gave specific batches to write state points
auto temp = get_node_array<int>(node_sp, "batches");
for (const auto& b : temp) {
statepoint_batch.insert(b);
}
} else {
// If neither were specified, write state point at last batch
statepoint_batch.insert(n_batches);
}
} else {
// If no <state_point> tag was present, by default write state point at
// last batch only
statepoint_batch.insert(n_batches);
}
// Check if the user has specified to write source points
if (check_for_node(root, "source_point")) {
// Get source_point node
xml_node node_sp = root.child("source_point");
// TODO: Read source point batches
// Determine batches at which to store source points
if (check_for_node(node_sp, "batches")) {
// User gave specific batches to write source points
auto temp = get_node_array<int>(node_sp, "batches");
for (const auto& b : temp) {
sourcepoint_batch.insert(b);
}
} else {
// If neither were specified, write source points with state points
sourcepoint_batch = statepoint_batch;
}
// Check if the user has specified to write binary source file
if (check_for_node(node_sp, "separate")) {
@ -609,10 +650,19 @@ void read_settings_xml()
// If no <source_point> tag was present, by default we keep source bank in
// statepoint file and write it out at statepoints intervals
source_separate = false;
// TODO: add defaults
sourcepoint_batch = statepoint_batch;
}
// TODO: Check source points are subset
// If source is not seperate and is to be written out in the statepoint file,
// make sure that the sourcepoint batch numbers are contained in the
// statepoint list
if (!source_separate) {
for (const auto& b : sourcepoint_batch) {
if (!contains(statepoint_batch, b)) {
fatal_error("Sourcepoint batches are not a subset of statepoint batches.");
}
}
}
// Check if the user has specified to not reduce tallies at the end of every
// batch
@ -778,6 +828,11 @@ extern "C" {
const char* openmc_path_particle_restart() {
return settings::path_particle_restart.c_str();
}
void free_memory_settings_c() {
settings::statepoint_batch.clear();
settings::sourcepoint_batch.clear();
}
}
} // namespace openmc

458
src/simulation.F90
View file

@ -7,392 +7,30 @@ module simulation
#endif
use bank_header, only: source_bank
use cmfd_execute, only: cmfd_init_batch, cmfd_tally_init, execute_cmfd
use cmfd_header, only: cmfd_on
use constants, only: ZERO
use eigenvalue, only: calculate_average_keff, calculate_generation_keff, &
k_sum
#ifdef _OPENMP
use eigenvalue, only: join_bank_from_threads
#endif
use error, only: fatal_error, write_message
use geometry_header, only: n_cells
use error, only: fatal_error
use material_header, only: n_materials, materials
use message_passing
use mgxs_interface, only: energy_bins, energy_bin_avg
use nuclide_header, only: micro_xs, n_nuclides
use output, only: header, print_columns, &
print_batch_keff, print_generation, print_runtime, &
print_results, write_tallies
use particle_header
use photon_header, only: micro_photon_xs, n_elements
use random_lcg, only: set_particle_seed
use settings
use simulation_header
use state_point, only: openmc_statepoint_write, write_source_point, load_state_point
use string, only: to_str
use tally, only: accumulate_tallies, setup_active_tallies, &
init_tally_routines
use tally_header
use tally_filter_header, only: filter_matches, n_filters, filter_match_pointer
use tally_derivative_header, only: tally_derivs
use timer_header
use trigger, only: check_triggers
use tracking, only: transport
implicit none
private
public :: openmc_next_batch
public :: openmc_simulation_init
public :: openmc_simulation_finalize
integer(C_INT), parameter :: STATUS_EXIT_NORMAL = 0
integer(C_INT), parameter :: STATUS_EXIT_MAX_BATCH = 1
integer(C_INT), parameter :: STATUS_EXIT_ON_TRIGGER = 2
interface
subroutine openmc_simulation_init_c() bind(C)
end subroutine
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
end function
end interface
contains
!===============================================================================
! OPENMC_NEXT_BATCH
!===============================================================================
function openmc_next_batch(status) result(err) bind(C)
integer(C_INT), intent(out), optional :: status
integer(C_INT) :: err
type(Particle) :: p
integer(8) :: i_work
err = 0
! Make sure simulation has been initialized
if (.not. simulation_initialized) then
err = E_ALLOCATE
call set_errmsg("Simulation has not been initialized yet.")
return
end if
call initialize_batch()
! =======================================================================
! LOOP OVER GENERATIONS
GENERATION_LOOP: do current_gen = 1, gen_per_batch
call initialize_generation()
! Start timer for transport
call time_transport % start()
! ====================================================================
! LOOP OVER PARTICLES
!$omp parallel do schedule(runtime) firstprivate(p) copyin(tally_derivs)
PARTICLE_LOOP: do i_work = 1, work
current_work = i_work
! grab source particle from bank
call initialize_history(p, current_work)
! transport particle
call transport(p)
end do PARTICLE_LOOP
!$omp end parallel do
! Accumulate time for transport
call time_transport % stop()
call finalize_generation()
end do GENERATION_LOOP
call finalize_batch()
! Check simulation ending criteria
if (present(status)) then
if (current_batch == n_max_batches) then
status = STATUS_EXIT_MAX_BATCH
elseif (satisfy_triggers) then
status = STATUS_EXIT_ON_TRIGGER
else
status = STATUS_EXIT_NORMAL
end if
end if
end function openmc_next_batch
!===============================================================================
! INITIALIZE_HISTORY
!===============================================================================
subroutine initialize_history(p, index_source)
type(Particle), intent(inout) :: p
integer(8), intent(in) :: index_source
integer(8) :: particle_seed ! unique index for particle
integer :: i
! set defaults
call particle_from_source(p, source_bank(index_source), run_CE, &
energy_bin_avg)
! set identifier for particle
p % id = work_index(rank) + index_source
! set random number seed
particle_seed = (total_gen + overall_generation() - 1)*n_particles + p % id
call set_particle_seed(particle_seed)
! set particle trace
trace = .false.
if (current_batch == trace_batch .and. current_gen == trace_gen .and. &
p % id == trace_particle) trace = .true.
! Set particle track.
p % write_track = .false.
if (write_all_tracks) then
p % write_track = .true.
else if (allocated(track_identifiers)) then
do i=1, size(track_identifiers(1,:))
if (current_batch == track_identifiers(1,i) .and. &
&current_gen == track_identifiers(2,i) .and. &
&p % id == track_identifiers(3,i)) then
p % write_track = .true.
exit
end if
end do
end if
end subroutine initialize_history
!===============================================================================
! INITIALIZE_BATCH
!===============================================================================
subroutine initialize_batch()
integer :: i
! Increment current batch
current_batch = current_batch + 1
if (run_mode == MODE_FIXEDSOURCE) then
call write_message("Simulating batch " // trim(to_str(current_batch)) &
// "...", 6)
end if
! Reset total starting particle weight used for normalizing tallies
total_weight = ZERO
if ((n_inactive > 0 .and. current_batch == 1) .or. &
(restart_run .and. restart_batch < n_inactive .and. current_batch == restart_batch + 1)) then
! Turn on inactive timer
call time_inactive % start()
elseif ((current_batch == n_inactive + 1) .or. &
(restart_run .and. restart_batch > n_inactive .and. current_batch == restart_batch + 1)) then
! Switch from inactive batch timer to active batch timer
call time_inactive % stop()
call time_active % start()
do i = 1, n_tallies
tallies(i) % obj % active = .true.
end do
end if
! check CMFD initialize batch
if (run_mode == MODE_EIGENVALUE) then
if (cmfd_run) call cmfd_init_batch()
end if
! Add user tallies to active tallies list
call setup_active_tallies()
end subroutine initialize_batch
!===============================================================================
! FINALIZE_GENERATION
!===============================================================================
subroutine finalize_generation()
interface
subroutine fill_source_bank_fixedsource() bind(C)
end subroutine
subroutine shannon_entropy() bind(C)
end subroutine
subroutine synchronize_bank() bind(C)
end subroutine
end interface
! Update global tallies with the omp private accumulation variables
!$omp parallel
!$omp critical
if (run_mode == MODE_EIGENVALUE) then
global_tallies(RESULT_VALUE, K_COLLISION) = &
global_tallies(RESULT_VALUE, K_COLLISION) + global_tally_collision
global_tallies(RESULT_VALUE, K_ABSORPTION) = &
global_tallies(RESULT_VALUE, K_ABSORPTION) + global_tally_absorption
global_tallies(RESULT_VALUE, K_TRACKLENGTH) = &
global_tallies(RESULT_VALUE, K_TRACKLENGTH) + global_tally_tracklength
end if
global_tallies(RESULT_VALUE, LEAKAGE) = &
global_tallies(RESULT_VALUE, LEAKAGE) + global_tally_leakage
!$omp end critical
! reset private tallies
if (run_mode == MODE_EIGENVALUE) then
global_tally_collision = ZERO
global_tally_absorption = ZERO
global_tally_tracklength = ZERO
end if
global_tally_leakage = ZERO
!$omp end parallel
if (run_mode == MODE_EIGENVALUE) then
#ifdef _OPENMP
! Join the fission bank from each thread into one global fission bank
call join_bank_from_threads()
#endif
! Distribute fission bank across processors evenly
call synchronize_bank()
! Calculate shannon entropy
if (entropy_on) call shannon_entropy()
! Collect results and statistics
call calculate_generation_keff()
call calculate_average_keff()
! Write generation output
if (master .and. verbosity >= 7) then
if (current_gen /= gen_per_batch) then
call print_generation()
end if
end if
elseif (run_mode == MODE_FIXEDSOURCE) then
! For fixed-source mode, we need to sample the external source
call fill_source_bank_fixedsource()
end if
end subroutine finalize_generation
!===============================================================================
! FINALIZE_BATCH handles synchronization and accumulation of tallies,
! calculation of Shannon entropy, getting single-batch estimate of keff, and
! turning on tallies when appropriate
!===============================================================================
subroutine finalize_batch()
integer(C_INT) :: err
character(MAX_FILE_LEN) :: filename
interface
subroutine broadcast_triggers() bind(C)
end subroutine broadcast_triggers
end interface
! Reduce tallies onto master process and accumulate
call time_tallies % start()
call accumulate_tallies()
call time_tallies % stop()
! Reset global tally results
if (current_batch <= n_inactive) then
global_tallies(:,:) = ZERO
n_realizations = 0
end if
if (run_mode == MODE_EIGENVALUE) then
! Perform CMFD calculation if on
if (cmfd_on) call execute_cmfd()
! Write batch output
if (master .and. verbosity >= 7) call print_batch_keff()
end if
! Check_triggers
if (master) call check_triggers()
#ifdef OPENMC_MPI
call broadcast_triggers()
#endif
if (satisfy_triggers .or. &
(trigger_on .and. current_batch == n_max_batches)) then
call statepoint_batch % add(current_batch)
end if
! Write out state point if it's been specified for this batch
if (statepoint_batch % contains(current_batch)) then
if (sourcepoint_batch % contains(current_batch) .and. source_write &
.and. .not. source_separate) then
err = openmc_statepoint_write(write_source=.true._C_BOOL)
else
err = openmc_statepoint_write(write_source=.false._C_BOOL)
end if
end if
! Write out a separate source point if it's been specified for this batch
if (sourcepoint_batch % contains(current_batch) .and. source_write &
.and. source_separate) call write_source_point()
! Write a continously-overwritten source point if requested.
if (source_latest) then
filename = trim(path_output) // 'source' // '.h5'
call write_source_point(filename)
end if
end subroutine finalize_batch
!===============================================================================
! INITIALIZE_SIMULATION
!===============================================================================
function openmc_simulation_init() result(err) bind(C)
integer(C_INT) :: err
subroutine simulation_init_f() bind(C)
integer :: i
err = 0
! Skip if simulation has already been initialized
if (simulation_initialized) return
! Call initialization on C++ side
call openmc_simulation_init_c()
! Set up tally procedure pointers
call init_tally_routines()
! Allocate source bank, and for eigenvalue simulations also allocate the
! fission bank
call allocate_banks()
! Allocate tally results arrays if they're not allocated yet
call configure_tallies()
! Activate the CMFD tallies
call cmfd_tally_init()
! Set up material nuclide index mapping
do i = 1, n_materials
call materials(i) % init_nuclide_index()
@ -410,66 +48,16 @@ contains
end do
!$omp end parallel
! 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 entropy_clear()
need_depletion_rx = .false.
! If this is a restart run, load the state point data and binary source
! file
if (restart_run) then
call load_state_point()
call write_message("Resuming simulation...", 6)
else
call initialize_source()
end if
! Display header
if (master) then
if (run_mode == MODE_FIXEDSOURCE) then
call header("FIXED SOURCE TRANSPORT SIMULATION", 3)
elseif (run_mode == MODE_EIGENVALUE) then
call header("K EIGENVALUE SIMULATION", 3)
if (verbosity >= 7) call print_columns()
end if
end if
! Set flag indicating initialization is done
simulation_initialized = .true.
end function openmc_simulation_init
end subroutine
!===============================================================================
! FINALIZE_SIMULATION calculates tally statistics, writes tallies, and displays
! execution time and results
!===============================================================================
function openmc_simulation_finalize() result(err) bind(C)
integer(C_INT) :: err
subroutine simulation_finalize_f() bind(C)
integer :: i ! loop index
interface
subroutine openmc_simulation_finalize_c() bind(C)
end subroutine openmc_simulation_finalize_c
subroutine print_overlap_check() bind(C)
end subroutine print_overlap_check
subroutine broadcast_results() bind(C)
end subroutine broadcast_results
end interface
err = 0
! Skip if simulation was never run
if (.not. simulation_initialized) return
! Stop active batch timer and start finalization timer
call time_active % stop()
call time_finalize % start()
integer :: i ! loop index
! Free up simulation-specific memory
do i = 1, n_materials
@ -479,45 +67,13 @@ contains
deallocate(micro_xs, micro_photon_xs, filter_matches)
!$omp end parallel
! Increment total number of generations
total_gen = total_gen + current_batch*gen_per_batch
#ifdef OPENMC_MPI
call broadcast_results()
#endif
! Write tally results to tallies.out
if (output_tallies .and. master) call write_tallies()
! Deactivate all tallies
if (allocated(tallies)) then
do i = 1, n_tallies
tallies(i) % obj % active = .false.
end do
end if
call openmc_simulation_finalize_c()
! Stop timers and show timing statistics
call time_finalize%stop()
call time_total%stop()
if (master) then
if (verbosity >= 6) call print_runtime()
if (verbosity >= 4) call print_results()
end if
if (check_overlaps) call print_overlap_check()
! Reset flags
need_depletion_rx = .false.
simulation_initialized = .false.
end function openmc_simulation_finalize
end subroutine
!===============================================================================
! ALLOCATE_BANKS allocates memory for the fission and source banks
!===============================================================================
subroutine allocate_banks()
subroutine allocate_banks() bind(C)
integer :: alloc_err ! allocation error code

415
src/simulation.cpp
View file

@ -1,13 +1,54 @@
#include "openmc/simulation.h"
#include "openmc/capi.h"
#include "openmc/container_util.h"
#include "openmc/eigenvalue.h"
#include "openmc/error.h"
#include "openmc/message_passing.h"
#include "openmc/output.h"
#include "openmc/particle.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/source.h"
#include "openmc/state_point.h"
#include "openmc/timer.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/tally.h"
#include <algorithm>
#include <string>
namespace openmc {
// data/functions from Fortran side
extern "C" bool cmfd_on;
extern "C" void accumulate_tallies();
extern "C" void allocate_banks();
extern "C" void check_triggers();
extern "C" void cmfd_init_batch();
extern "C" void cmfd_tally_init();
extern "C" void configure_tallies();
extern "C" void execute_cmfd();
extern "C" void init_tally_routines();
extern "C" void join_bank_from_threads();
extern "C" void load_state_point();
extern "C" void print_batch_keff();
extern "C" void print_columns();
extern "C" void print_generation();
extern "C" void print_results();
extern "C" void print_runtime();
extern "C" void setup_active_tallies();
extern "C" void simulation_init_f();
extern "C" void simulation_finalize_f();
extern "C" void transport(Particle* p);
extern "C" void write_tallies();
} // namespace openmc
//==============================================================================
// C API functions
//==============================================================================
// OPENMC_RUN encompasses all the main logic where iterations are performed
// over the batches, generations, and histories in a fixed source or k-eigenvalue
@ -27,6 +68,176 @@ int openmc_run()
return err;
}
int openmc_simulation_init()
{
using namespace openmc;
// Skip if simulation has already been initialized
if (simulation::initialized) return 0;
// Determine how much work each process should do
calculate_work();
// Allocate array for matching filter bins
#pragma omp parallel
{
filter_matches.resize(n_filters);
}
// Set up tally procedure pointers
init_tally_routines();
// Allocate source bank, and for eigenvalue simulations also allocate the
// fission bank
allocate_banks();
// Allocate tally results arrays if they're not allocated yet
configure_tallies();
// Activate the CMFD tallies
cmfd_tally_init();
// Call Fortran initialization
simulation_init_f();
// Reset global variables -- this is done before loading state point (as that
// will potentially populate k_generation and entropy)
simulation::current_batch = 0;
simulation::k_generation.clear();
entropy.clear();
simulation::need_depletion_rx = false;
// If this is a restart run, load the state point data and binary source
// file
if (settings::restart_run) {
load_state_point();
write_message("Resuming simulation...", 6);
} else {
initialize_source();
}
// Display header
if (mpi::master) {
if (settings::run_mode == RUN_MODE_FIXEDSOURCE) {
header("FIXED SOURCE TRANSPORT SIMULATION", 3);
} else if (settings::run_mode == RUN_MODE_EIGENVALUE) {
header("K EIGENVALUE SIMULATION", 3);
if (settings::verbosity >= 7) print_columns();
}
}
// Set flag indicating initialization is done
simulation::initialized = true;
return 0;
}
int openmc_simulation_finalize()
{
using namespace openmc;
// Skip if simulation was never run
if (!simulation::initialized) return 0;
// Stop active batch timer and start finalization timer
time_active.stop();
time_finalize.start();
#pragma omp parallel
{
filter_matches.clear();
}
// Deallocate Fortran variables, set tallies to inactive
simulation_finalize_f();
// Increment total number of generations
simulation::total_gen += simulation::current_batch*settings::gen_per_batch;
#ifdef OPENMC_MPI
broadcast_results();
#endif
// Write tally results to tallies.out
if (settings::output_tallies && mpi::master) write_tallies();
// Deactivate all tallies
for (int i = 1; i <= n_tallies; ++i) {
openmc_tally_set_active(i, false);
}
// Stop timers and show timing statistics
time_finalize.stop();
time_total.stop();
if (mpi::master) {
if (settings::verbosity >= 6) print_runtime();
if (settings::verbosity >= 4) print_results();
}
if (settings::check_overlaps) print_overlap_check();
// Reset flags
simulation::need_depletion_rx = false;
simulation::initialized = false;
return 0;
}
int openmc_next_batch(int* status)
{
using namespace openmc;
using openmc::simulation::current_gen;
// Make sure simulation has been initialized
if (!simulation::initialized) {
set_errmsg("Simulation has not been initialized yet.");
return OPENMC_E_ALLOCATE;
}
initialize_batch();
// =======================================================================
// LOOP OVER GENERATIONS
for (current_gen = 1; current_gen <= settings::gen_per_batch; ++current_gen) {
initialize_generation();
// Start timer for transport
time_transport.start();
// ====================================================================
// LOOP OVER PARTICLES
#pragma omp parallel for schedule(runtime)
for (int64_t i_work = 1; i_work <= simulation::work; ++i_work) {
simulation::current_work = i_work;
// grab source particle from bank
Particle p;
initialize_history(&p, simulation::current_work);
// transport particle
transport(&p);
}
// Accumulate time for transport
time_transport.stop();
finalize_generation();
}
finalize_batch();
// Check simulation ending criteria
if (status) {
if (simulation::current_batch == settings::n_max_batches) {
*status = STATUS_EXIT_MAX_BATCH;
} else if (simulation::satisfy_triggers) {
*status = STATUS_EXIT_ON_TRIGGER;
} else {
*status = STATUS_EXIT_NORMAL;
}
}
return 0;
}
namespace openmc {
//==============================================================================
@ -38,6 +249,7 @@ namespace simulation {
int current_batch;
int current_gen;
int64_t current_work;
bool initialized {false};
double keff {1.0};
double keff_std;
double k_col_abs {0.0};
@ -48,7 +260,6 @@ int n_lost_particles {0};
bool need_depletion_rx {false};
int restart_batch;
bool satisfy_triggers {false};
bool simulation_initialized {false};
int total_gen {0};
int64_t work;
@ -68,15 +279,104 @@ int thread_id; //!< ID of a given thread
// Non-member functions
//==============================================================================
void openmc_simulation_init_c()
void initialize_batch()
{
// Determine how much work each process should do
calculate_work();
// Increment current batch
++simulation::current_batch;
// Allocate array for matching filter bins
#pragma omp parallel
{
filter_matches.resize(n_filters);
if (settings::run_mode == RUN_MODE_FIXEDSOURCE) {
int b = simulation::current_batch;
write_message("Simulating batch " + std::to_string(b), 6);
}
// Reset total starting particle weight used for normalizing tallies
total_weight = 0.0;
// Determine if this batch is the first inactive or active batch.
bool first_inactive = false;
bool first_active = false;
if (!settings::restart_run) {
first_inactive = settings::n_inactive > 0 && simulation::current_batch == 1;
first_active = simulation::current_batch == settings::n_inactive + 1;
} else if (simulation::current_batch == simulation::restart_batch + 1){
first_inactive = simulation::restart_batch < settings::n_inactive;
first_active = !first_inactive;
}
// Manage active/inactive timers and activate tallies if necessary.
if (first_inactive) {
time_inactive.start();
} else if (first_active) {
time_inactive.stop();
time_active.start();
for (int i = 1; i <= n_tallies; ++i) {
// TODO: change one-based index
openmc_tally_set_active(i, true);
}
}
// check CMFD initialize batch
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
if (settings::cmfd_run) cmfd_init_batch();
}
// Add user tallies to active tallies list
setup_active_tallies();
}
void finalize_batch()
{
// Reduce tallies onto master process and accumulate
time_tallies.start();
accumulate_tallies();
time_tallies.stop();
// Reset global tally results
if (simulation::current_batch <= settings::n_inactive) {
auto gt = global_tallies();
std::fill(gt.begin(), gt.end(), 0.0);
n_realizations = 0;
}
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
// Perform CMFD calculation if on
if (cmfd_on) execute_cmfd();
// Write batch output
if (mpi::master && settings::verbosity >= 7) print_batch_keff();
}
// Check_triggers
if (mpi::master) check_triggers();
#ifdef OPENMC_MPI
MPI_Bcast(&simulation::satisfy_triggers, 1, MPI_C_BOOL, 0, mpi::intracomm);
#endif
if (simulation::satisfy_triggers || (settings::trigger_on &&
simulation::current_batch == settings::n_max_batches)) {
settings::statepoint_batch.insert(simulation::current_batch);
}
// Write out state point if it's been specified for this batch
if (contains(settings::statepoint_batch, simulation::current_batch)) {
if (contains(settings::sourcepoint_batch, simulation::current_batch)
&& settings::source_write && !settings::source_separate) {
bool b = true;
openmc_statepoint_write(nullptr, &b);
} else {
bool b = false;
openmc_statepoint_write(nullptr, &b);
}
}
// Write out a separate source point if it's been specified for this batch
if (contains(settings::sourcepoint_batch, simulation::current_batch)
&& settings::source_write && settings::source_separate) {
write_source_point(nullptr);
}
// Write a continously-overwritten source point if requested.
if (settings::source_latest) {
auto filename = settings::path_output + "source.h5";
write_source_point(filename.c_str());
}
}
@ -94,12 +394,99 @@ void initialize_generation()
}
}
extern "C" void
openmc_simulation_finalize_c()
void finalize_generation()
{
#pragma omp parallel
auto gt = global_tallies();
// Update global tallies with the omp private accumulation variables
#pragma omp parallel
{
filter_matches.clear();
#pragma omp critical(increment_global_tallies)
{
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
gt(K_COLLISION, RESULT_VALUE) += global_tally_collision;
gt(K_ABSORPTION, RESULT_VALUE) += global_tally_absorption;
gt(K_TRACKLENGTH, RESULT_VALUE) += global_tally_tracklength;
}
gt(LEAKAGE, RESULT_VALUE) += global_tally_leakage;
}
// reset threadprivate tallies
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
global_tally_collision = 0.0;
global_tally_absorption = 0.0;
global_tally_tracklength = 0.0;
}
global_tally_leakage = 0.0;
}
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
#ifdef _OPENMP
// Join the fission bank from each thread into one global fission bank
join_bank_from_threads();
#endif
// Distribute fission bank across processors evenly
synchronize_bank();
// Calculate shannon entropy
if (settings::entropy_on) shannon_entropy();
// Collect results and statistics
calculate_generation_keff();
calculate_average_keff();
// Write generation output
if (mpi::master && settings::verbosity >= 7) {
if (simulation::current_gen != settings::gen_per_batch) {
print_generation();
}
}
} else if (settings::run_mode == RUN_MODE_FIXEDSOURCE) {
// For fixed-source mode, we need to sample the external source
fill_source_bank_fixedsource();
}
}
void initialize_history(Particle* p, int64_t index_source)
{
// Get pointer to source bank
Bank* source_bank;
int64_t n;
openmc_source_bank(&source_bank, &n);
// set defaults
p->from_source(&source_bank[index_source - 1]);
// set identifier for particle
p->id = simulation::work_index[mpi::rank] + index_source;
// set random number seed
int64_t particle_seed = (simulation::total_gen + overall_generation() - 1)
* settings::n_particles + p->id;
set_particle_seed(particle_seed);
// set particle trace
simulation::trace = false;
if (simulation::current_batch == settings::trace_batch &&
simulation::current_gen == settings::trace_gen &&
p->id == settings::trace_particle) simulation::trace = true;
// Set particle track.
p->write_track = false;
if (settings::write_all_tracks) {
p->write_track = true;
} else if (settings::track_identifiers.size() > 0) {
for (const auto& t : settings::track_identifiers) {
if (simulation::current_batch == t[0] &&
simulation::current_gen == t[1] &&
p->id == t[2]) {
p->write_track = true;
break;
}
}
}
}
@ -165,10 +552,6 @@ void broadcast_results() {
simulation::k_abs_tra = temp[2];
}
void broadcast_triggers()
{
MPI_Bcast(&simulation::satisfy_triggers, 1, MPI_C_BOOL, 0, mpi::intracomm);
}
#endif
//==============================================================================

9
src/simulation_header.F90
View file

@ -23,7 +23,6 @@ module simulation_header
integer(C_INT), bind(C) :: current_batch ! current batch
integer(C_INT), bind(C) :: current_gen ! current generation within a batch
integer(C_INT), bind(C) :: total_gen ! total number of generations simulated
logical(C_BOOL), bind(C) :: simulation_initialized
logical(C_BOOL), bind(C) :: need_depletion_rx ! need to calculate depletion reaction rx?
! ============================================================================
@ -32,7 +31,6 @@ module simulation_header
logical(C_BOOL), bind(C) :: satisfy_triggers ! whether triggers are satisfied
integer(C_INT64_T), bind(C) :: work ! number of particles per processor
integer(C_INT64_T), bind(C) :: current_work ! index in source bank of current history simulated
! ============================================================================
! K-EIGENVALUE SIMULATION VARIABLES
@ -59,7 +57,7 @@ module simulation_header
logical(C_BOOL), bind(C) :: trace
!$omp threadprivate(trace, thread_id, current_work)
!$omp threadprivate(trace, thread_id)
interface
subroutine entropy_clear() bind(C)
@ -84,11 +82,6 @@ module simulation_header
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
function work_index(rank) result(i) bind(C)
import C_INT, C_INT64_T
integer(C_INT), value :: rank

6
src/source.cpp
View file

@ -352,11 +352,7 @@ extern "C" double total_source_strength()
return strength;
}
// Needed in fill_source_bank_fixedsource
extern "C" int overall_generation();
//! Fill source bank at end of generation for fixed source simulations
extern "C" void fill_source_bank_fixedsource()
void fill_source_bank_fixedsource()
{
if (settings::path_source.empty()) {
// Get pointer to source bank

71
src/state_point.F90
View file

@ -16,7 +16,6 @@ module state_point
use bank_header, only: Bank
use cmfd_header
use constants
use eigenvalue, only: openmc_get_keff, k_sum
use endf, only: reaction_name
use error, only: fatal_error, warning, write_message
use hdf5_interface
@ -391,19 +390,19 @@ contains
! Write out the runtime metrics.
runtime_group = create_group(file_id, "runtime")
call write_dataset(runtime_group, "total initialization", &
time_initialize % get_value())
time_initialize_elapsed())
call write_dataset(runtime_group, "reading cross sections", &
time_read_xs % get_value())
call write_dataset(runtime_group, "simulation", &
time_inactive % get_value() + time_active % get_value())
time_inactive_elapsed() + time_active_elapsed())
call write_dataset(runtime_group, "transport", &
time_transport % get_value())
time_transport_elapsed())
if (run_mode == MODE_EIGENVALUE) then
call write_dataset(runtime_group, "inactive batches", &
time_inactive % get_value())
time_inactive_elapsed())
end if
call write_dataset(runtime_group, "active batches", &
time_active % get_value())
time_active_elapsed())
if (run_mode == MODE_EIGENVALUE) then
call write_dataset(runtime_group, "synchronizing fission bank", &
time_bank_elapsed())
@ -413,7 +412,7 @@ contains
time_bank_sendrecv_elapsed())
end if
call write_dataset(runtime_group, "accumulating tallies", &
time_tallies % get_value())
time_tallies_elapsed())
if (cmfd_run) then
call write_dataset(runtime_group, "CMFD", time_cmfd % get_value())
call write_dataset(runtime_group, "CMFD building matrices", &
@ -421,7 +420,7 @@ contains
call write_dataset(runtime_group, "CMFD solving matrices", &
time_cmfdsolve % get_value())
end if
call write_dataset(runtime_group, "total", time_total % get_value())
call write_dataset(runtime_group, "total", time_total_elapsed())
call close_group(runtime_group)
call file_close(file_id)
@ -443,52 +442,13 @@ contains
end if
end function openmc_statepoint_write
!===============================================================================
! WRITE_SOURCE_POINT
!===============================================================================
subroutine write_source_point(filename)
character(MAX_FILE_LEN), intent(in), optional :: filename
logical :: parallel
integer(HID_T) :: file_id
character(MAX_FILE_LEN) :: filename_
! When using parallel HDF5, the file is written to collectively by all
! processes. With MPI-only, the file is opened and written by the master
! (note that the call to write_source_bank is by all processes since slave
! processes need to send source bank data to the master.
#ifdef PHDF5
parallel = .true.
#else
parallel = .false.
#endif
if (present(filename)) then
filename_ = filename
else
filename_ = trim(path_output) // 'source.' // &
& zero_padded(current_batch, count_digits(n_max_batches))
filename_ = trim(filename_) // '.h5'
end if
if (master .or. parallel) then
file_id = file_open(filename_, 'w', parallel=.true.)
call write_attribute(file_id, "filetype", 'source')
end if
call write_source_bank(file_id, source_bank)
if (master .or. parallel) call file_close(file_id)
end subroutine write_source_point
!===============================================================================
! LOAD_STATE_POINT
!===============================================================================
subroutine load_state_point()
subroutine load_state_point() bind(C)
integer :: i
integer :: n
integer :: int_array(3)
integer, allocatable :: array(:)
integer(C_INT64_T) :: seed
@ -504,6 +464,9 @@ contains
import HID_T
integer(HID_T), value :: group
end subroutine
subroutine restart_set_keff() bind(C)
end subroutine
end interface
! Write message
@ -612,17 +575,7 @@ contains
! Set k_sum, keff, and current_batch based on whether restart file is part
! 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(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
n = k_generation_size()
keff = k_generation(n)
end if
call restart_set_keff()
current_batch = restart_batch
! Check to make sure source bank is present

57
src/state_point.cpp
View file

@ -1,6 +1,7 @@
#include "openmc/state_point.h"
#include <algorithm>
#include <iomanip> // for setfill, setw
#include <string>
#include <vector>
@ -9,6 +10,7 @@
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/eigenvalue.h"
#include "openmc/error.h"
#include "openmc/hdf5_interface.h"
#include "openmc/message_passing.h"
@ -18,7 +20,6 @@
namespace openmc {
hid_t h5banktype() {
// Create type for array of 3 reals
hsize_t dims[] {3};
@ -36,6 +37,46 @@ hid_t h5banktype() {
return banktype;
}
void
write_source_point(const char* filename)
{
// When using parallel HDF5, the file is written to collectively by all
// processes. With MPI-only, the file is opened and written by the master
// (note that the call to write_source_bank is by all processes since slave
// processes need to send source bank data to the master.
#ifdef PHDF5
bool parallel = true;
#else
bool parallel = false;
#endif
std::string filename_;
if (filename) {
filename_ = filename;
} else {
// Determine width for zero padding
int w = std::to_string(settings::n_max_batches).size();
std::stringstream s;
s << settings::path_output << "source." << std::setfill('0')
<< std::setw(w) << simulation::current_batch << ".h5";
filename_ = s.str();
}
hid_t file_id;
if (mpi::master || parallel) {
file_id = file_open(filename_, 'w', true);
write_attribute(file_id, "filetype", "source");
}
// Get pointer to source bank and write to file
Bank* source_bank;
int64_t n;
openmc_source_bank(&source_bank, &n);
write_source_bank(file_id, source_bank);
if (mpi::master || parallel) file_close(file_id);
}
void
write_source_bank(hid_t group_id, Bank* source_bank)
@ -275,4 +316,18 @@ void write_tally_results_nr(hid_t file_id)
}
}
void restart_set_keff()
{
if (simulation::restart_batch > settings::n_inactive) {
for (int i = settings::n_inactive; i < simulation::restart_batch; ++i) {
k_sum[0] += simulation::k_generation[i];
k_sum[1] += std::pow(simulation::k_generation[i], 2);
}
int n = settings::gen_per_batch*n_realizations;
simulation::keff = k_sum[0] / n;
} else {
simulation::keff = simulation::k_generation.back();
}
}
} // namespace openmc

6
src/tallies/tally.F90
View file

@ -52,7 +52,7 @@ contains
! with the CE and MG modes.
!===============================================================================
subroutine init_tally_routines()
subroutine init_tally_routines() bind(C)
if (run_CE) then
score_general => score_general_ce
score_analog_tally => score_analog_tally_ce
@ -3765,7 +3765,7 @@ contains
! within the batch to a new random variable
!===============================================================================
subroutine accumulate_tallies()
subroutine accumulate_tallies() bind(C)
integer :: i
real(C_DOUBLE) :: k_col ! Copy of batch collision estimate of keff
@ -3826,7 +3826,7 @@ contains
! SETUP_ACTIVE_TALLIES
!===============================================================================
subroutine setup_active_tallies()
subroutine setup_active_tallies() bind(C)
integer :: i ! loop counter

9
src/tallies/tally.cpp
View file

@ -13,6 +13,15 @@
namespace openmc {
//==============================================================================
// Global variable definitions
//==============================================================================
double global_tally_absorption;
double global_tally_collision;
double global_tally_tracklength;
double global_tally_leakage;
//==============================================================================
// Non-member functions
//==============================================================================

10
src/tallies/tally_header.F90
View file

@ -136,10 +136,10 @@ module tally_header
! global tally, it can cause a higher cache miss rate due to
! invalidation. Thus, we use threadprivate variables to accumulate global
! tallies and then reduce at the end of a generation.
real(C_DOUBLE), public :: global_tally_collision = ZERO
real(C_DOUBLE), public, bind(C) :: global_tally_absorption = ZERO
real(C_DOUBLE), public :: global_tally_tracklength = ZERO
real(C_DOUBLE), public :: global_tally_leakage = ZERO
real(C_DOUBLE), public, bind(C) :: global_tally_collision
real(C_DOUBLE), public, bind(C) :: global_tally_absorption
real(C_DOUBLE), public, bind(C) :: global_tally_tracklength
real(C_DOUBLE), public, bind(C) :: global_tally_leakage
!$omp threadprivate(global_tally_collision, global_tally_absorption, &
!$omp& global_tally_tracklength, global_tally_leakage)
@ -392,7 +392,7 @@ contains
! tallies.xml file.
!===============================================================================
subroutine configure_tallies()
subroutine configure_tallies() bind(C)
integer :: i

2
src/tallies/trigger.F90
View file

@ -29,7 +29,7 @@ contains
! and predicts the number of remainining batches to convergence.
!===============================================================================
subroutine check_triggers()
subroutine check_triggers() bind(C)
implicit none

27
src/timer.cpp
View file

@ -38,23 +38,48 @@ double Timer::elapsed()
// Global variables
//==============================================================================
Timer time_active;
Timer time_bank;
Timer time_bank_sample;
Timer time_bank_sendrecv;
Timer time_finalize;
Timer time_inactive;
Timer time_initialize;
Timer time_tallies;
Timer time_total;
Timer time_transport;
//==============================================================================
// Fortran compatibility
//==============================================================================
extern "C" double time_active_elapsed() { return time_active.elapsed(); }
extern "C" double time_bank_elapsed() { return time_bank.elapsed(); }
extern "C" double time_bank_sample_elapsed() { return time_bank_sample.elapsed(); }
extern "C" double time_bank_sendrecv_elapsed() { return time_bank_sendrecv.elapsed(); }
extern "C" double time_finalize_elapsed() { return time_finalize.elapsed(); }
extern "C" double time_inactive_elapsed() { return time_inactive.elapsed(); }
extern "C" double time_initialize_elapsed() { return time_initialize.elapsed(); }
extern "C" double time_tallies_elapsed() { return time_tallies.elapsed(); }
extern "C" double time_total_elapsed() { return time_total.elapsed(); }
extern "C" double time_transport_elapsed() { return time_transport.elapsed(); }
extern "C" void reset_timers()
//==============================================================================
// Non-member functions
//==============================================================================
void reset_timers()
{
time_active.reset();
time_bank.reset();
time_bank_sample.reset();
time_bank_sendrecv.reset();
time_finalize.reset();
time_inactive.reset();
time_initialize.reset();
time_tallies.reset();
time_total.reset();
time_transport.reset();
}
} // namespace openmc

44
src/timer_header.F90
View file

@ -7,6 +7,10 @@ module timer_header
implicit none
interface
function time_active_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_bank_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
@ -19,6 +23,30 @@ module timer_header
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_finalize_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_inactive_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_initialize_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_tallies_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_total_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_transport_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
subroutine reset_timers() bind(C)
end subroutine
end interface
@ -44,15 +72,7 @@ module timer_header
! ============================================================================
! TIMING VARIABLES
type(Timer) :: time_total ! timer for total run
type(Timer) :: time_initialize ! timer for initialization
type(Timer) :: time_read_xs ! timer for reading cross sections
type(Timer) :: time_unionize ! timer for material xs-energy grid union
type(Timer) :: time_tallies ! timer for accumulate tallies
type(Timer) :: time_inactive ! timer for inactive batches
type(Timer) :: time_active ! timer for active batches
type(Timer) :: time_transport ! timer for transport only
type(Timer) :: time_finalize ! timer for finalization
contains
@ -117,4 +137,12 @@ contains
self % elapsed = ZERO
end subroutine timer_reset
!===============================================================================
! RESET_TIMERS resets timers on the Fortran side
!===============================================================================
subroutine reset_timers_f() bind(C)
call time_read_xs % reset()
end subroutine
end module timer_header

5
src/tracking.F90
View file

@ -48,7 +48,7 @@ contains
! TRANSPORT encompasses the main logic for moving a particle through geometry.
!===============================================================================
subroutine transport(p)
subroutine transport(p) bind(C)
type(Particle), intent(inout) :: p
@ -287,8 +287,7 @@ contains
! Check for secondary particles if this particle is dead
if (.not. p % alive) then
if (p % n_secondary > 0) then
call particle_from_source(p, p % secondary_bank(p % n_secondary), &
run_CE, energy_bin_avg)
call particle_from_source(p, p % secondary_bank(p % n_secondary))
p % n_secondary = p % n_secondary - 1
n_event = 0