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