diff --git a/CMakeLists.txt b/CMakeLists.txt index dfd9221ce..913b9a00c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -251,9 +251,26 @@ elseif(CMAKE_C_COMPILER_ID MATCHES Clang) endif() +list(APPEND cxxflags -std=c++11 -O2) +if(debug) + list(REMOVE_ITEM cxxflags -O2) + list(APPEND cxxflags -g -O0) +endif() +if(profile) + list(APPEND cxxflags -pg) +endif() +if(optimize) + list(REMOVE_ITEM cxxflags -O2) + list(APPEND cxxflags -O3) +endif() +if(openmp) + list(APPEND cxxflags -fopenmp) +endif() + # Show flags being used message(STATUS "Fortran flags: ${f90flags}") message(STATUS "C flags: ${cflags}") +message(STATUS "C++ flags: ${cxxflags}") message(STATUS "Linker flags: ${ldflags}") #=============================================================================== @@ -312,7 +329,7 @@ endif() add_library(faddeeva STATIC src/faddeeva/Faddeeva.c) #=============================================================================== -# Build OpenMC executable +# List source files. Define the libopenmc and the OpenMC executable #=============================================================================== set(program "openmc") @@ -418,19 +435,37 @@ set(LIBOPENMC_FORTRAN_SRC src/tallies/trigger.F90 src/tallies/trigger_header.F90 ) -add_library(libopenmc SHARED ${LIBOPENMC_FORTRAN_SRC}) +set(LIBOPENMC_CXX_SRC + src/random_lcg.h + src/random_lcg.cpp) +add_library(libopenmc SHARED ${LIBOPENMC_FORTRAN_SRC} ${LIBOPENMC_CXX_SRC}) set_target_properties(libopenmc PROPERTIES OUTPUT_NAME openmc) add_executable(${program} src/main.F90) + +#=============================================================================== +# Add compiler/linker flags +#=============================================================================== + set_property(TARGET ${program} libopenmc pugixml_fortran PROPERTY LINKER_LANGUAGE Fortran) target_include_directories(libopenmc PUBLIC ${HDF5_INCLUDE_DIRS}) -# set compile flags via target_compile_options +# The executable and the faddeeva package use only one language. They can be +# set via target_compile_options which accepts a list. target_compile_options(${program} PUBLIC ${f90flags}) -target_compile_options(libopenmc PUBLIC ${f90flags}) target_compile_options(faddeeva PRIVATE ${cflags}) +# The libopenmc library has both F90 and C++ so the compile flags must be set +# file-by-file via set_source_file_properties. The compile flags must first be +# converted from lists to strings. +string(REPLACE ";" " " f90flags "${f90flags}") +string(REPLACE ";" " " cxxflags "${cxxflags}") +set_source_files_properties(${LIBOPENMC_FORTRAN_SRC} PROPERTIES COMPILE_FLAGS + ${f90flags}) +set_source_files_properties(${LIBOPENMC_CXX_SRC} PROPERTIES COMPILE_FLAGS + ${cxxflags}) + # Add HDF5 library directories to link line with -L foreach(LIBDIR ${HDF5_LIBRARY_DIRS}) list(APPEND ldflags "-L${LIBDIR}") @@ -438,7 +473,8 @@ endforeach() # target_link_libraries treats any arguments starting with - but not -l as # linker flags. Thus, we can pass both linker flags and libraries together. -target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} pugixml_fortran faddeeva) +target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} pugixml_fortran + faddeeva) target_link_libraries(${program} ${ldflags} libopenmc) #=============================================================================== diff --git a/src/constants.F90 b/src/constants.F90 index 73aca25d2..5bc9eb88e 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -1,5 +1,7 @@ module constants + use, intrinsic :: ISO_C_BINDING + implicit none ! ============================================================================ @@ -413,12 +415,13 @@ module constants ! ============================================================================ ! RANDOM NUMBER STREAM CONSTANTS - integer, parameter :: N_STREAMS = 5 - integer, parameter :: STREAM_TRACKING = 1 - integer, parameter :: STREAM_TALLIES = 2 - integer, parameter :: STREAM_SOURCE = 3 - integer, parameter :: STREAM_URR_PTABLE = 4 - integer, parameter :: STREAM_VOLUME = 5 + integer(C_INT), bind(C, name='N_STREAMS') :: N_STREAMS + integer(C_INT), bind(C, name='STREAM_TRACKING') :: STREAM_TRACKING + integer(C_INT), bind(C, name='STREAM_TALLIES') :: STREAM_TALLIES + integer(C_INT), bind(C, name='STREAM_SOURCE') :: STREAM_SOURCE + integer(C_INT), bind(C, name='STREAM_URR_PTABLE') :: STREAM_URR_PTABLE + integer(C_INT), bind(C, name='STREAM_VOLUME') :: STREAM_VOLUME + integer(C_INT64_T), parameter :: DEFAULT_SEED = 1_8 ! ============================================================================ ! MISCELLANEOUS CONSTANTS diff --git a/src/initialize.F90 b/src/initialize.F90 index 117381328..21d568c72 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -20,7 +20,7 @@ module initialize use message_passing use mgxs_data, only: read_mgxs, create_macro_xs use output, only: print_version, write_message, print_usage - use random_lcg, only: openmc_set_seed, seed + use random_lcg, only: openmc_set_seed use settings #ifdef _OPENMP use simulation_header, only: n_threads @@ -84,7 +84,7 @@ contains ! Initialize random number generator -- if the user specifies a seed, it ! will be re-initialized later - err = openmc_set_seed(seed) + err = openmc_set_seed(DEFAULT_SEED) ! Read XML input files call read_input_xml() diff --git a/src/random_lcg.F90 b/src/random_lcg.F90 index 685557f07..839bd729a 100644 --- a/src/random_lcg.F90 +++ b/src/random_lcg.F90 @@ -2,199 +2,47 @@ module random_lcg use, intrinsic :: ISO_C_BINDING - use constants - implicit none - private - save + integer(C_INT64_T), bind(C) :: seed - ! Starting seed - integer(C_INT64_T), public, bind(C) :: seed = 1_8 + interface + function prn() result(pseudo_rn) bind(C) + use ISO_C_BINDING + implicit none + real(C_DOUBLE) :: pseudo_rn + end function prn - ! LCG parameters - integer(C_INT64_T), parameter :: prn_mult = 2806196910506780709_8 ! multiplication factor, g - integer(C_INT64_T), parameter :: prn_add = 1_8 ! additive factor, c - integer, parameter :: prn_bits = 63 ! number of bits, M - integer(C_INT64_T), parameter :: prn_mod = ibset(0_8, prn_bits) ! 2^M - integer(C_INT64_T), parameter :: prn_mask = not(prn_mod) ! 2^M - 1 - integer(C_INT64_T), parameter :: prn_stride = 152917_8 ! stride between particles - real(C_DOUBLE), parameter :: prn_norm = 2._8**(-prn_bits) ! 2^(-M) + function future_prn(n) result(pseudo_rn) bind(C) + use ISO_C_BINDING + implicit none + integer(C_INT64_T), value :: n + real(C_DOUBLE) :: pseudo_rn + end function future_prn - ! Current PRNG state - integer(C_INT64_T) :: prn_seed(N_STREAMS) ! current seed - integer :: stream ! current RNG stream -!$omp threadprivate(prn_seed, stream) + subroutine set_particle_seed(id) bind(C) + use ISO_C_BINDING + implicit none + integer(C_INT64_T), value :: id + end subroutine set_particle_seed - public :: prn - public :: future_prn - public :: set_particle_seed - public :: advance_prn_seed - public :: prn_set_stream - public :: openmc_set_seed + subroutine advance_prn_seed(n) bind(C) + use ISO_C_BINDING + implicit none + integer(C_INT64_T), value :: n + end subroutine advance_prn_seed -contains - -!=============================================================================== -! PRN generates a pseudo-random number using a linear congruential generator -!=============================================================================== - - function prn() result(pseudo_rn) - - real(C_DOUBLE) :: pseudo_rn - - ! This algorithm uses bit-masking to find the next integer(C_INT64_T) value - ! to be used to calculate the random number - - prn_seed(stream) = iand(prn_mult*prn_seed(stream) + prn_add, prn_mask) - - ! Once the integer is calculated, we just need to divide by 2**m, - ! represented here as multiplying by a pre-calculated factor - - pseudo_rn = prn_seed(stream) * prn_norm - - end function prn - -!=============================================================================== -! FUTURE_PRN generates a pseudo-random number which is 'n' times ahead from the -! current seed. -!=============================================================================== - - function future_prn(n) result(pseudo_rn) - - integer(C_INT64_T), intent(in) :: n ! number of prns to skip - - real(C_DOUBLE) :: pseudo_rn - - pseudo_rn = future_seed(n, prn_seed(stream)) * prn_norm - - end function future_prn - -!=============================================================================== -! SET_PARTICLE_SEED sets the seed to a unique value based on the ID of the -! particle -!=============================================================================== - - subroutine set_particle_seed(id) - - integer(C_INT64_T), intent(in) :: id - - integer :: i - - do i = 1, N_STREAMS - prn_seed(i) = future_seed(id*prn_stride, seed + i - 1) - end do - - end subroutine set_particle_seed - -!=============================================================================== -! ADVANCE_PRN_SEED advances the random number seed 'n' times from the current -! seed. -!=============================================================================== - - subroutine advance_prn_seed(n) - - integer(C_INT64_T), intent(in) :: n ! number of seeds to skip - - prn_seed(stream) = future_seed(n, prn_seed(stream)) - - end subroutine advance_prn_seed - -!=============================================================================== -! FUTURE_SEED advances the random number seed 'skip' times. This is usually -! used to skip a fixed number of random numbers (the stride) so that a given -! particle always has the same starting seed regardless of how many processors -! are used -!=============================================================================== - - function future_seed(n, seed) result(new_seed) - - integer(C_INT64_T), intent(in) :: n ! number of seeds to skip - integer(C_INT64_T), intent(in) :: seed ! original seed - integer(C_INT64_T) :: new_seed ! new seed - - integer(C_INT64_T) :: nskip ! positive number of seeds to skip - integer(C_INT64_T) :: g ! original multiplicative constant - integer(C_INT64_T) :: c ! original additive constnat - integer(C_INT64_T) :: g_new ! new effective multiplicative constant - integer(C_INT64_T) :: c_new ! new effective additive constant - - ! In cases where we want to skip backwards, we add the period of the random - ! number generator until the number of PRNs to skip is positive since - ! skipping ahead that much is the same as skipping backwards by the original - ! amount - - nskip = n - do while (nskip < 0_8) - nskip = nskip + prn_mod - end do - - ! Make sure nskip is less than 2^M - nskip = iand(nskip, prn_mask) - - ! The algorithm here to determine the parameters used to skip ahead is - ! described in F. Brown, "Random Number Generation with Arbitrary Stride," - ! Trans. Am. Nucl. Soc. (Nov. 1994). This algorithm is able to skip ahead in - ! O(log2(N)) operations instead of O(N). Basically, it computes parameters G - ! and C which can then be used to find x_N = G*x_0 + C mod 2^M. - - ! Initialize constants - g = prn_mult - c = prn_add - g_new = 1 - c_new = 0 - BIT_LOOP: do while (nskip > 0_8) - ! Check if least significant bit is 1 - if (btest(nskip,0)) then - g_new = iand(g_new*g, prn_mask) - c_new = iand(c_new*g + c, prn_mask) - endif - c = iand((g+1)*c, prn_mask) - g = iand(g*g, prn_mask) - - ! Move bits right, dropping least significant bit - nskip = ishft(nskip, -1) - end do BIT_LOOP - - ! With G and C, we can now find the new seed - new_seed = iand(g_new*seed + c_new, prn_mask) - - end function future_seed - -!=============================================================================== -! PRN_SET_STREAM changes the random number stream. If random numbers are needed -! in routines not used directly for tracking (e.g. physics), this allows the -! numbers to be generated without affecting reproducibility of the physics. -!=============================================================================== - - subroutine prn_set_stream(i) - - integer, intent(in) :: i - - stream = i - - end subroutine prn_set_stream - -!=============================================================================== -! C API FUNCTIONS -!=============================================================================== - - function openmc_set_seed(new_seed) result(err) bind(C) - ! Saves the starting seed and sets up the PRNG thread state - integer(C_INT64_T), value, intent(in) :: new_seed - integer(C_INT) :: err - - integer :: i - - err = 0 - seed = new_seed -!$omp parallel - do i = 1, N_STREAMS - prn_seed(i) = seed + i - 1 - end do - stream = STREAM_TRACKING -!$omp end parallel - - end function openmc_set_seed + subroutine prn_set_stream(n) bind(C) + use ISO_C_BINDING + implicit none + integer(C_INT), value :: n + end subroutine prn_set_stream + function openmc_set_seed(new_seed) result(err) bind(C) + use ISO_C_BINDING + implicit none + integer(C_INT64_T), value :: new_seed + integer(C_INT) :: err + end function openmc_set_seed + end interface end module random_lcg diff --git a/src/random_lcg.cpp b/src/random_lcg.cpp new file mode 100644 index 000000000..047bda8e5 --- /dev/null +++ b/src/random_lcg.cpp @@ -0,0 +1,145 @@ +#include "random_lcg.h" +#include + + +// Constants +extern "C" const int N_STREAMS {5}; +extern "C" const int STREAM_TRACKING {0}; +extern "C" const int STREAM_TALLIES {1}; +extern "C" const int STREAM_SOURCE {2}; +extern "C" const int STREAM_URR_PTABLE {3}; +extern "C" const int STREAM_VOLUME {4}; + +// Starting seed +int64_t seed {1}; + +// LCG parameters +constexpr uint64_t prn_mult {2806196910506780709LL}; // multiplication + // factor, g +constexpr uint64_t prn_add {1}; // additive factor, c +constexpr uint64_t prn_mod {0x8000000000000000}; // 2^63 +constexpr uint64_t prn_mask {0x7fffffffffffffff}; // 2^63 - 1 +constexpr uint64_t prn_stride {152917LL}; // stride between + // particles +constexpr double prn_norm {1.0 / prn_mod}; // 2^-63 + +// Current PRNG state +uint64_t prn_seed[N_STREAMS]; // current seed +int stream; // current RNG stream +#pragma omp threadprivate(prn_seed, stream) + + +//============================================================================== +// PRN +//============================================================================== + +extern "C" double +prn() +{ + // This algorithm uses bit-masking to find the next integer(8) value to be + // used to calculate the random number. + prn_seed[stream] = (prn_mult*prn_seed[stream] + prn_add) & prn_mask; + + // Once the integer is calculated, we just need to divide by 2**m, + // represented here as multiplying by a pre-calculated factor + return prn_seed[stream] * prn_norm; +} + +//============================================================================== +// FUTURE_PRN +//============================================================================== + +extern "C" double +future_prn(int64_t n) +{ + return future_seed(static_cast(n), prn_seed[stream]) * prn_norm; +} + +//============================================================================== +// SET_PARTICLE_SEED +//============================================================================== + +extern "C" void +set_particle_seed(int64_t id) +{ + for (int i = 0; i < N_STREAMS; i++) { + prn_seed[i] = future_seed(static_cast(id) * prn_stride, seed + i); + } +} + +//============================================================================== +// ADVANCE_PRN_SEED +//============================================================================== + +extern "C" void +advance_prn_seed(int64_t n) +{ + prn_seed[stream] = future_seed(static_cast(n), prn_seed[stream]); +} + +//============================================================================== +// FUTURE_SEED +//============================================================================== + +uint64_t +future_seed(uint64_t n, uint64_t seed) +{ + // Make sure nskip is less than 2^M. + n &= prn_mask; + + // The algorithm here to determine the parameters used to skip ahead is + // described in F. Brown, "Random Number Generation with Arbitrary Stride," + // Trans. Am. Nucl. Soc. (Nov. 1994). This algorithm is able to skip ahead in + // O(log2(N)) operations instead of O(N). Basically, it computes parameters G + // and C which can then be used to find x_N = G*x_0 + C mod 2^M. + + // Initialize constants + uint64_t g {prn_mult}; + uint64_t c {prn_add}; + uint64_t g_new {1}; + uint64_t c_new {0}; + + while (n > 0) { + // Check if the least significant bit is 1. + if (n & 1) { + g_new *= g; + c_new = c_new * g + c; + } + c *= (g + 1); + g *= g; + + // Move bits right, dropping least significant bit. + n >>= 1; + } + + // With G and C, we can now find the new seed. + return (g_new * seed + c_new) & prn_mask; +} + +//============================================================================== +// PRN_SET_STREAM +//============================================================================== + +extern "C" void +prn_set_stream(int i) +{ + stream = i; // Shift by one to move from Fortran to C indexing. +} + +//============================================================================== +// API FUNCTIONS +//============================================================================== + +extern "C" int +openmc_set_seed(int64_t new_seed) +{ + seed = new_seed; +#pragma omp parallel + { + for (int i = 0; i < N_STREAMS; i++) { + prn_seed[i] = seed + i; + } + prn_set_stream(STREAM_TRACKING); + } + return 0; +} diff --git a/src/random_lcg.h b/src/random_lcg.h new file mode 100644 index 000000000..04191254f --- /dev/null +++ b/src/random_lcg.h @@ -0,0 +1,85 @@ +#ifndef RANDOM_LCG_H +#define RANDOM_LCG_H + +#include + +//============================================================================== +// Module constants. +//============================================================================== + +extern "C" const int N_STREAMS; +extern "C" const int STREAM_TRACKING; +extern "C" const int STREAM_TALLIES; +extern "C" const int STREAM_SOURCE; +extern "C" const int STREAM_URR_PTABLE; +extern "C" const int STREAM_VOLUME; + +//============================================================================== +//! Generate a pseudo-random number using a linear congruential generator. +//! @return A random number between 0 and 1 +//============================================================================== + +extern "C" double prn(); + +//============================================================================== +//! Generate a random number which is 'n' times ahead from the current seed. +//! +//! The result of this function will be the same as the result from calling +//! `prn()` 'n' times. +//! @param n The number of RNG seeds to skip ahead by +//! @return A random number between 0 and 1 +//============================================================================== + +extern "C" double future_prn(int64_t n); + +//============================================================================== +//! Set the RNG seed to a unique value based on the ID of the particle. +//! @param id The particle ID +//============================================================================== + +extern "C" void set_particle_seed(int64_t id); + +//============================================================================== +//! Advance the random number seed 'n' times from the current seed. +//! @param n The number of RNG seeds to skip ahead by +//============================================================================== + +extern "C" void advance_prn_seed(int64_t n); + +//============================================================================== +//! Advance a random number seed 'n' times. +//! +//! This is usually used to skip a fixed number of random numbers (the stride) +//! so that a given particle always has the same starting seed regardless of +//! how many processors are used. +//! @param n The number of RNG seeds to skip ahead by +//! @param seed The starting to seed to advance from +//============================================================================== + +uint64_t future_seed(uint64_t n, uint64_t seed); + +//============================================================================== +//! Switch the RNG to a different stream of random numbers. +//! +//! If random numbers are needed in routines not used directly for tracking +//! (e.g. physics), this allows the numbers to be generated without affecting +//! reproducibility of the physics. +//! @param n The RNG stream to switch to. Use the constants such as +//! `STREAM_TRACKING` and `STREAM_TALLIES` for this argument. +//============================================================================== + +extern "C" void prn_set_stream(int n); + +//============================================================================== +// API FUNCTIONS +//============================================================================== + +//============================================================================== +//! Set OpenMC's master seed. +//! @param new_seed The master seed. All other seeds will be derived from this +//! one. +//============================================================================== + +extern "C" int openmc_set_seed(int64_t new_seed); + +#endif // RANDOM_LCG_H