From cc52bfc435f6d00ea6939ceb5ef6931840005780 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 1 Dec 2017 16:45:28 -0500 Subject: [PATCH 1/8] Translate the RNG to C++ --- CMakeLists.txt | 46 +++++++++- src/random_lcg.C | 153 +++++++++++++++++++++++++++++++ src/random_lcg.F90 | 223 ++++++++------------------------------------- src/random_lcg.h | 56 ++++++++++++ 4 files changed, 286 insertions(+), 192 deletions(-) create mode 100644 src/random_lcg.C create mode 100644 src/random_lcg.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 0b0bef14e..a8b801eb7 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 "CXX flags: ${cxxflags}") message(STATUS "Linker flags: ${ldflags}") #=============================================================================== @@ -292,7 +309,7 @@ set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE) add_library(faddeeva STATIC src/faddeeva/Faddeeva.c) #=============================================================================== -# Build OpenMC executable +# List source files. Define the libopenmc and the OpenMC executable #=============================================================================== set(program "openmc") @@ -398,19 +415,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.C) +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}") @@ -418,7 +453,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/random_lcg.C b/src/random_lcg.C new file mode 100644 index 000000000..c78c40748 --- /dev/null +++ b/src/random_lcg.C @@ -0,0 +1,153 @@ +#include "random_lcg.h" +#include + + +// Starting seed +int64_t seed = 1; + +// LCG parameters +const int64_t prn_mult = 2806196910506780709LL; // multiplication factor, g +const int64_t prn_add = 1; // additive factor, c +const int64_t prn_mod = -0x8000000000000000; // -2^63 +const int64_t prn_mask = 0x7fffffffffffffff; // 2^63 - 1 +const int64_t prn_stride = 152917LL; // stride between particles +const double prn_norm = pow(2, -63); // 2^-63 + +// Module constants +const int N_STREAMS = 5; +const int STREAM_TRACKING = 0; +const int STREAM_TALLIES = 1; +const int STREAM_SOURCE = 2; +const int STREAM_URR_PTABLE = 3; +const int STREAM_VOLUME = 4; + +// Current PRNG state +int64_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 + double pseudo_rn = prn_seed[stream] * prn_norm; + return pseudo_rn; +} + +//============================================================================== +// FUTURE_PRN +//============================================================================== + +extern "C" double +future_prn(int64_t n) +{ + double pseudo_rn = future_seed(n, prn_seed[stream]) * prn_norm; + return pseudo_rn; +} + +//============================================================================== +// 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(id * prn_stride, seed + i); + } +} + +//============================================================================== +// ADVANCE_PRN_SEED +//============================================================================== + +extern "C" void +advance_prn_seed(int64_t n) +{ + prn_seed[stream] = future_seed(n, prn_seed[stream]); +} + +//============================================================================== +// FUTURE_SEED +//============================================================================== + +extern "C" int64_t +future_seed(int64_t n, int64_t seed) +{ + // 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. + + int64_t nskip = n; + while (nskip < 0) nskip += prn_mod; + + // Make sure nskip is less than 2^M. + 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 + int64_t g = prn_mult; + int64_t c = prn_add; + int64_t g_new = 1; + int64_t c_new = 0; + + while (nskip > 0) { + // Check if the least significant bit is 1. + if (nskip & 1) { + g_new = (g_new * g) & prn_mask; + c_new = (c_new * g + c) & prn_mask; + } + c = ((g + 1) * c) & prn_mask; + g = (g * g) & prn_mask; + + // Move bits right, dropping least significant bit. + nskip >>= 1; + } + + // With G and C, we can now find the new seed. + int64_t new_seed = (g_new * seed + c_new) & prn_mask; + return new_seed; +} + +//============================================================================== +// PRN_SET_STREAM +//============================================================================== + +extern "C" void +prn_set_stream(int i) +{ + stream = i - 1; +} + +//============================================================================== +// 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; + } + stream = STREAM_TRACKING; +#pragma end omp parallel + return 0; +} diff --git a/src/random_lcg.F90 b/src/random_lcg.F90 index 685557f07..b76341cbd 100644 --- a/src/random_lcg.F90 +++ b/src/random_lcg.F90 @@ -2,199 +2,48 @@ 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, name='prn') + 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, name='future_prn') + 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, name='set_particle_seed') + 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, name='advance_prn_seed') + 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, name='prn_set_stream') + 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, name='openmc_set_seed') + 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.h b/src/random_lcg.h new file mode 100644 index 000000000..90d1c253e --- /dev/null +++ b/src/random_lcg.h @@ -0,0 +1,56 @@ +#ifndef RANDOM_LCG_H +#define RANDOM_LCG_H + +#include + +//============================================================================== +// PRN generates a pseudo-random number using a linear congruential generator. +//============================================================================== + +extern "C" double prn(); + +//============================================================================== +// FUTURE_PRN generates a pseudo-random number which is 'n' times ahead from the +// current seed. +//============================================================================== + +extern "C" double future_prn(int64_t n); + +//============================================================================== +// SET_PARTICLE_SEED sets the seed to a unique value based on the ID of the +// particle. +//============================================================================== + +extern "C" void set_particle_seed(int64_t id); + +//============================================================================== +// ADVANCE_PRN_SEED advances the random number seed 'n' times from the current +// seed. +//============================================================================== + +extern "C" void advance_prn_seed(int64_t n); + +//============================================================================== +// 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. +//============================================================================== + +extern "C" int64_t future_seed(int64_t n, int64_t 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. +//============================================================================== + +extern "C" void prn_set_stream(int n); + +//============================================================================== +// API FUNCTIONS +//============================================================================== + +extern "C" int openmc_set_seed(int64_t new_seed); + +#endif // RANDOM_LCG_H From 6811a197325a094db06a051860c85daf438b2435 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 1 Dec 2017 17:32:27 -0500 Subject: [PATCH 2/8] Use unsigned integers in the RNG --- src/random_lcg.C | 49 +++++++++++++++++++++++------------------------- src/random_lcg.h | 10 +++++----- 2 files changed, 28 insertions(+), 31 deletions(-) diff --git a/src/random_lcg.C b/src/random_lcg.C index c78c40748..8383401f9 100644 --- a/src/random_lcg.C +++ b/src/random_lcg.C @@ -6,12 +6,12 @@ int64_t seed = 1; // LCG parameters -const int64_t prn_mult = 2806196910506780709LL; // multiplication factor, g -const int64_t prn_add = 1; // additive factor, c -const int64_t prn_mod = -0x8000000000000000; // -2^63 -const int64_t prn_mask = 0x7fffffffffffffff; // 2^63 - 1 -const int64_t prn_stride = 152917LL; // stride between particles -const double prn_norm = pow(2, -63); // 2^-63 +const uint64_t prn_mult = 2806196910506780709LL; // multiplication factor, g +const uint64_t prn_add = 1; // additive factor, c +const uint64_t prn_mod = 0x8000000000000000; // 2^63 +const uint64_t prn_mask = 0x7fffffffffffffff; // 2^63 - 1 +const uint64_t prn_stride = 152917LL; // stride between particles +const double prn_norm = pow(2, -63); // 2^-63 // Module constants const int N_STREAMS = 5; @@ -22,8 +22,8 @@ const int STREAM_URR_PTABLE = 3; const int STREAM_VOLUME = 4; // Current PRNG state -int64_t prn_seed[N_STREAMS]; // current seed -int stream; // current RNG stream +uint64_t prn_seed[N_STREAMS]; // current seed +int stream; // current RNG stream #pragma omp threadprivate(prn_seed, stream) @@ -40,8 +40,7 @@ prn() // Once the integer is calculated, we just need to divide by 2**m, // represented here as multiplying by a pre-calculated factor - double pseudo_rn = prn_seed[stream] * prn_norm; - return pseudo_rn; + return prn_seed[stream] * prn_norm; } //============================================================================== @@ -49,10 +48,9 @@ prn() //============================================================================== extern "C" double -future_prn(int64_t n) +future_prn(uint64_t n) { - double pseudo_rn = future_seed(n, prn_seed[stream]) * prn_norm; - return pseudo_rn; + return future_seed(n, prn_seed[stream]) * prn_norm; } //============================================================================== @@ -60,7 +58,7 @@ future_prn(int64_t n) //============================================================================== extern "C" void -set_particle_seed(int64_t id) +set_particle_seed(uint64_t id) { for (int i = 0; i < N_STREAMS; i++) { prn_seed[i] = future_seed(id * prn_stride, seed + i); @@ -72,7 +70,7 @@ set_particle_seed(int64_t id) //============================================================================== extern "C" void -advance_prn_seed(int64_t n) +advance_prn_seed(uint64_t n) { prn_seed[stream] = future_seed(n, prn_seed[stream]); } @@ -81,16 +79,16 @@ advance_prn_seed(int64_t n) // FUTURE_SEED //============================================================================== -extern "C" int64_t -future_seed(int64_t n, int64_t seed) +extern "C" uint64_t +future_seed(uint64_t n, uint64_t seed) { // 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. - int64_t nskip = n; - while (nskip < 0) nskip += prn_mod; + uint64_t nskip = n; + while (nskip > prn_mod) nskip += prn_mod; // Make sure nskip is less than 2^M. nskip &= prn_mask; @@ -102,10 +100,10 @@ future_seed(int64_t n, int64_t seed) // and C which can then be used to find x_N = G*x_0 + C mod 2^M. // Initialize constants - int64_t g = prn_mult; - int64_t c = prn_add; - int64_t g_new = 1; - int64_t c_new = 0; + uint64_t g = prn_mult; + uint64_t c = prn_add; + uint64_t g_new = 1; + uint64_t c_new = 0; while (nskip > 0) { // Check if the least significant bit is 1. @@ -121,8 +119,7 @@ future_seed(int64_t n, int64_t seed) } // With G and C, we can now find the new seed. - int64_t new_seed = (g_new * seed + c_new) & prn_mask; - return new_seed; + return (g_new * seed + c_new) & prn_mask; } //============================================================================== @@ -140,7 +137,7 @@ prn_set_stream(int i) //============================================================================== extern "C" int -openmc_set_seed(int64_t new_seed) +openmc_set_seed(uint64_t new_seed) { seed = new_seed; #pragma omp parallel diff --git a/src/random_lcg.h b/src/random_lcg.h index 90d1c253e..889195f10 100644 --- a/src/random_lcg.h +++ b/src/random_lcg.h @@ -14,21 +14,21 @@ extern "C" double prn(); // current seed. //============================================================================== -extern "C" double future_prn(int64_t n); +extern "C" double future_prn(uint64_t n); //============================================================================== // SET_PARTICLE_SEED sets the seed to a unique value based on the ID of the // particle. //============================================================================== -extern "C" void set_particle_seed(int64_t id); +extern "C" void set_particle_seed(uint64_t id); //============================================================================== // ADVANCE_PRN_SEED advances the random number seed 'n' times from the current // seed. //============================================================================== -extern "C" void advance_prn_seed(int64_t n); +extern "C" void advance_prn_seed(uint64_t n); //============================================================================== // FUTURE_SEED advances the random number seed 'skip' times. This is usually @@ -37,7 +37,7 @@ extern "C" void advance_prn_seed(int64_t n); // are used. //============================================================================== -extern "C" int64_t future_seed(int64_t n, int64_t seed); +extern "C" uint64_t future_seed(uint64_t n, uint64_t seed); //============================================================================== // PRN_SET_STREAM changes the random number stream. If random numbers are needed @@ -51,6 +51,6 @@ extern "C" void prn_set_stream(int n); // API FUNCTIONS //============================================================================== -extern "C" int openmc_set_seed(int64_t new_seed); +extern "C" int openmc_set_seed(uint64_t new_seed); #endif // RANDOM_LCG_H From 1cfa7c89568450f110e570dbf57b2266bd2d1a37 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 4 Dec 2017 15:19:49 -0500 Subject: [PATCH 3/8] Simplify RNG --- src/random_lcg.C | 28 ++++++++++------------------ 1 file changed, 10 insertions(+), 18 deletions(-) diff --git a/src/random_lcg.C b/src/random_lcg.C index 8383401f9..ef49cd13b 100644 --- a/src/random_lcg.C +++ b/src/random_lcg.C @@ -11,10 +11,10 @@ const uint64_t prn_add = 1; // additive factor, c const uint64_t prn_mod = 0x8000000000000000; // 2^63 const uint64_t prn_mask = 0x7fffffffffffffff; // 2^63 - 1 const uint64_t prn_stride = 152917LL; // stride between particles -const double prn_norm = pow(2, -63); // 2^-63 +const double prn_norm = pow(2, -63); // 2^-63 // Module constants -const int N_STREAMS = 5; +const int N_STREAMS = 5; const int STREAM_TRACKING = 0; const int STREAM_TALLIES = 1; const int STREAM_SOURCE = 2; @@ -82,16 +82,8 @@ advance_prn_seed(uint64_t n) extern "C" uint64_t future_seed(uint64_t n, uint64_t seed) { - // 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. - - uint64_t nskip = n; - while (nskip > prn_mod) nskip += prn_mod; - // Make sure nskip is less than 2^M. - nskip &= prn_mask; + 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," @@ -105,17 +97,17 @@ future_seed(uint64_t n, uint64_t seed) uint64_t g_new = 1; uint64_t c_new = 0; - while (nskip > 0) { + while (n > 0) { // Check if the least significant bit is 1. - if (nskip & 1) { - g_new = (g_new * g) & prn_mask; - c_new = (c_new * g + c) & prn_mask; + if (n & 1) { + g_new = g_new * g; + c_new = c_new * g + c; } - c = ((g + 1) * c) & prn_mask; - g = (g * g) & prn_mask; + c = (g + 1) * c; + g = g * g; // Move bits right, dropping least significant bit. - nskip >>= 1; + n >>= 1; } // With G and C, we can now find the new seed. From 168f202e3ecf48cdd15b541dc396b24465832986 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 11 Dec 2017 15:31:31 -0500 Subject: [PATCH 4/8] Address #939 feedback --- CMakeLists.txt | 2 +- src/constants.F90 | 3 +++ src/initialize.F90 | 4 ++-- src/random_lcg.C | 28 ++++++++++++++-------------- src/random_lcg.h | 8 ++++---- 5 files changed, 24 insertions(+), 21 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index a8b801eb7..41ca8d8ee 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -270,7 +270,7 @@ endif() # Show flags being used message(STATUS "Fortran flags: ${f90flags}") message(STATUS "C flags: ${cflags}") -message(STATUS "CXX flags: ${cxxflags}") +message(STATUS "C++ flags: ${cxxflags}") message(STATUS "Linker flags: ${ldflags}") #=============================================================================== diff --git a/src/constants.F90 b/src/constants.F90 index 73aca25d2..2466dbb14 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -1,5 +1,7 @@ module constants + use, intrinsic :: ISO_C_BINDING + implicit none ! ============================================================================ @@ -419,6 +421,7 @@ module constants integer, parameter :: STREAM_SOURCE = 3 integer, parameter :: STREAM_URR_PTABLE = 4 integer, parameter :: STREAM_VOLUME = 5 + 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.C b/src/random_lcg.C index ef49cd13b..5bdf13d70 100644 --- a/src/random_lcg.C +++ b/src/random_lcg.C @@ -15,11 +15,11 @@ const double prn_norm = pow(2, -63); // 2^-63 // Module constants const int N_STREAMS = 5; -const int STREAM_TRACKING = 0; -const int STREAM_TALLIES = 1; -const int STREAM_SOURCE = 2; -const int STREAM_URR_PTABLE = 3; -const int STREAM_VOLUME = 4; +const int STREAM_TRACKING = 1; +const int STREAM_TALLIES = 2; +const int STREAM_SOURCE = 3; +const int STREAM_URR_PTABLE = 4; +const int STREAM_VOLUME = 5; // Current PRNG state uint64_t prn_seed[N_STREAMS]; // current seed @@ -48,9 +48,9 @@ prn() //============================================================================== extern "C" double -future_prn(uint64_t n) +future_prn(int64_t n) { - return future_seed(n, prn_seed[stream]) * prn_norm; + return future_seed(static_cast(n), prn_seed[stream]) * prn_norm; } //============================================================================== @@ -58,10 +58,10 @@ future_prn(uint64_t n) //============================================================================== extern "C" void -set_particle_seed(uint64_t id) +set_particle_seed(int64_t id) { for (int i = 0; i < N_STREAMS; i++) { - prn_seed[i] = future_seed(id * prn_stride, seed + i); + prn_seed[i] = future_seed(static_cast(id) * prn_stride, seed + i); } } @@ -70,16 +70,16 @@ set_particle_seed(uint64_t id) //============================================================================== extern "C" void -advance_prn_seed(uint64_t n) +advance_prn_seed(int64_t n) { - prn_seed[stream] = future_seed(n, prn_seed[stream]); + prn_seed[stream] = future_seed(static_cast(n), prn_seed[stream]); } //============================================================================== // FUTURE_SEED //============================================================================== -extern "C" uint64_t +uint64_t future_seed(uint64_t n, uint64_t seed) { // Make sure nskip is less than 2^M. @@ -121,7 +121,7 @@ future_seed(uint64_t n, uint64_t seed) extern "C" void prn_set_stream(int i) { - stream = i - 1; + stream = i - 1; // Shift by one to move from Fortran to C indexing. } //============================================================================== @@ -136,7 +136,7 @@ openmc_set_seed(uint64_t new_seed) for (int i = 0; i < N_STREAMS; i++) { prn_seed[i] = seed + i; } - stream = STREAM_TRACKING; + prn_set_stream(STREAM_TRACKING); #pragma end omp parallel return 0; } diff --git a/src/random_lcg.h b/src/random_lcg.h index 889195f10..4c7ec0140 100644 --- a/src/random_lcg.h +++ b/src/random_lcg.h @@ -14,21 +14,21 @@ extern "C" double prn(); // current seed. //============================================================================== -extern "C" double future_prn(uint64_t n); +extern "C" double future_prn(int64_t n); //============================================================================== // SET_PARTICLE_SEED sets the seed to a unique value based on the ID of the // particle. //============================================================================== -extern "C" void set_particle_seed(uint64_t id); +extern "C" void set_particle_seed(int64_t id); //============================================================================== // ADVANCE_PRN_SEED advances the random number seed 'n' times from the current // seed. //============================================================================== -extern "C" void advance_prn_seed(uint64_t n); +extern "C" void advance_prn_seed(int64_t n); //============================================================================== // FUTURE_SEED advances the random number seed 'skip' times. This is usually @@ -37,7 +37,7 @@ extern "C" void advance_prn_seed(uint64_t n); // are used. //============================================================================== -extern "C" uint64_t future_seed(uint64_t n, uint64_t seed); +uint64_t future_seed(uint64_t n, uint64_t seed); //============================================================================== // PRN_SET_STREAM changes the random number stream. If random numbers are needed From ed5ff626f4ff2c35eb524c89198be7f6eb73e213 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 12 Dec 2017 15:16:54 -0500 Subject: [PATCH 5/8] Link Fortran RNG constants to C++ definitions --- src/constants.F90 | 12 ++++++------ src/random_lcg.C | 12 ++---------- src/random_lcg.F90 | 13 ++++++------- src/random_lcg.h | 13 ++++++++++++- 4 files changed, 26 insertions(+), 24 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index 2466dbb14..5bc9eb88e 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -415,12 +415,12 @@ 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 ! ============================================================================ diff --git a/src/random_lcg.C b/src/random_lcg.C index 5bdf13d70..42450fd37 100644 --- a/src/random_lcg.C +++ b/src/random_lcg.C @@ -13,14 +13,6 @@ const uint64_t prn_mask = 0x7fffffffffffffff; // 2^63 - 1 const uint64_t prn_stride = 152917LL; // stride between particles const double prn_norm = pow(2, -63); // 2^-63 -// Module constants -const int N_STREAMS = 5; -const int STREAM_TRACKING = 1; -const int STREAM_TALLIES = 2; -const int STREAM_SOURCE = 3; -const int STREAM_URR_PTABLE = 4; -const int STREAM_VOLUME = 5; - // Current PRNG state uint64_t prn_seed[N_STREAMS]; // current seed int stream; // current RNG stream @@ -121,7 +113,7 @@ future_seed(uint64_t n, uint64_t seed) extern "C" void prn_set_stream(int i) { - stream = i - 1; // Shift by one to move from Fortran to C indexing. + stream = i; // Shift by one to move from Fortran to C indexing. } //============================================================================== @@ -129,7 +121,7 @@ prn_set_stream(int i) //============================================================================== extern "C" int -openmc_set_seed(uint64_t new_seed) +openmc_set_seed(int64_t new_seed) { seed = new_seed; #pragma omp parallel diff --git a/src/random_lcg.F90 b/src/random_lcg.F90 index b76341cbd..839bd729a 100644 --- a/src/random_lcg.F90 +++ b/src/random_lcg.F90 @@ -7,39 +7,38 @@ module random_lcg integer(C_INT64_T), bind(C) :: seed interface - function prn() result(pseudo_rn) bind(C, name='prn') + function prn() result(pseudo_rn) bind(C) use ISO_C_BINDING implicit none real(C_DOUBLE) :: pseudo_rn end function prn - function future_prn(n) result(pseudo_rn) bind(C, name='future_prn') + 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 - subroutine set_particle_seed(id) bind(C, name='set_particle_seed') + subroutine set_particle_seed(id) bind(C) use ISO_C_BINDING implicit none integer(C_INT64_T), value :: id end subroutine set_particle_seed - subroutine advance_prn_seed(n) bind(C, name='advance_prn_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 - subroutine prn_set_stream(n) bind(C, name='prn_set_stream') + 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, name='openmc_set_seed') + function openmc_set_seed(new_seed) result(err) bind(C) use ISO_C_BINDING implicit none integer(C_INT64_T), value :: new_seed diff --git a/src/random_lcg.h b/src/random_lcg.h index 4c7ec0140..ef8ed8a44 100644 --- a/src/random_lcg.h +++ b/src/random_lcg.h @@ -3,6 +3,17 @@ #include +//============================================================================== +// Module 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; + //============================================================================== // PRN generates a pseudo-random number using a linear congruential generator. //============================================================================== @@ -51,6 +62,6 @@ extern "C" void prn_set_stream(int n); // API FUNCTIONS //============================================================================== -extern "C" int openmc_set_seed(uint64_t new_seed); +extern "C" int openmc_set_seed(int64_t new_seed); #endif // RANDOM_LCG_H From ec302b52d64befdd53d4a0ac17daa16c935e660b Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 12 Dec 2017 17:57:53 -0500 Subject: [PATCH 6/8] Fix multiple definition error in random_lcg.h --- src/random_lcg.C | 38 ++++++++++++++++++++++++-------------- src/random_lcg.h | 12 ++++++------ 2 files changed, 30 insertions(+), 20 deletions(-) diff --git a/src/random_lcg.C b/src/random_lcg.C index 42450fd37..ea79be67f 100644 --- a/src/random_lcg.C +++ b/src/random_lcg.C @@ -2,16 +2,26 @@ #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; +int64_t seed {1}; // LCG parameters -const uint64_t prn_mult = 2806196910506780709LL; // multiplication factor, g -const uint64_t prn_add = 1; // additive factor, c -const uint64_t prn_mod = 0x8000000000000000; // 2^63 -const uint64_t prn_mask = 0x7fffffffffffffff; // 2^63 - 1 -const uint64_t prn_stride = 152917LL; // stride between particles -const double prn_norm = pow(2, -63); // 2^-63 +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 {pow(2, -63)}; // 2^-63 // Current PRNG state uint64_t prn_seed[N_STREAMS]; // current seed @@ -84,19 +94,19 @@ future_seed(uint64_t n, uint64_t seed) // 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; + 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_new * g; + g_new *= g; c_new = c_new * g + c; } - c = (g + 1) * c; - g = g * g; + c *= (g + 1); + g *= g; // Move bits right, dropping least significant bit. n >>= 1; diff --git a/src/random_lcg.h b/src/random_lcg.h index ef8ed8a44..ca8dfd826 100644 --- a/src/random_lcg.h +++ b/src/random_lcg.h @@ -7,12 +7,12 @@ // Module 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; +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; //============================================================================== // PRN generates a pseudo-random number using a linear congruential generator. From f6d6b926b15d70b1ff7477f4bc8d7ea46504ff05 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 13 Dec 2017 18:49:12 -0500 Subject: [PATCH 7/8] Address #939 comments --- CMakeLists.txt | 2 +- src/{random_lcg.C => random_lcg.cpp} | 11 ++++++----- 2 files changed, 7 insertions(+), 6 deletions(-) rename src/{random_lcg.C => random_lcg.cpp} (95%) diff --git a/CMakeLists.txt b/CMakeLists.txt index 41ca8d8ee..4edc89ebc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -417,7 +417,7 @@ set(LIBOPENMC_FORTRAN_SRC ) set(LIBOPENMC_CXX_SRC src/random_lcg.h - src/random_lcg.C) + 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) diff --git a/src/random_lcg.C b/src/random_lcg.cpp similarity index 95% rename from src/random_lcg.C rename to src/random_lcg.cpp index ea79be67f..047bda8e5 100644 --- a/src/random_lcg.C +++ b/src/random_lcg.cpp @@ -21,7 +21,7 @@ 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 {pow(2, -63)}; // 2^-63 +constexpr double prn_norm {1.0 / prn_mod}; // 2^-63 // Current PRNG state uint64_t prn_seed[N_STREAMS]; // current seed @@ -135,10 +135,11 @@ 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; + { + for (int i = 0; i < N_STREAMS; i++) { + prn_seed[i] = seed + i; + } + prn_set_stream(STREAM_TRACKING); } - prn_set_stream(STREAM_TRACKING); -#pragma end omp parallel return 0; } From 4b6402104e0cbccc9cfdd23aeec78e85f5e8d595 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 14 Dec 2017 16:19:26 -0500 Subject: [PATCH 8/8] Switch C++ RNG to Doxygen-style comments --- src/random_lcg.h | 46 ++++++++++++++++++++++++++++++++-------------- 1 file changed, 32 insertions(+), 14 deletions(-) diff --git a/src/random_lcg.h b/src/random_lcg.h index ca8dfd826..04191254f 100644 --- a/src/random_lcg.h +++ b/src/random_lcg.h @@ -15,45 +15,57 @@ extern "C" const int STREAM_URR_PTABLE; extern "C" const int STREAM_VOLUME; //============================================================================== -// PRN generates a pseudo-random number using a linear congruential generator. +//! Generate a pseudo-random number using a linear congruential generator. +//! @return A random number between 0 and 1 //============================================================================== extern "C" double prn(); //============================================================================== -// FUTURE_PRN generates a pseudo-random number which is 'n' times ahead from the -// current seed. +//! 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_PARTICLE_SEED sets the seed to a unique value based on the ID of the -// particle. +//! 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_PRN_SEED advances the random number seed 'n' times from the current -// seed. +//! 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); //============================================================================== -// 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. +//! 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); //============================================================================== -// 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. +//! 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); @@ -62,6 +74,12 @@ 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