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

Translate the RNG to C++

Browse source
This commit is contained in:
Sterling Harper 2017-12-01 16:45:28 -05:00
parent 75c7d41022
commit cc52bfc435
4 changed files with 286 additions and 192 deletions
Download patch file Download diff file Expand all files Collapse all files

46
CMakeLists.txt
View file

@ -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)
#===============================================================================

153
src/random_lcg.C Normal file
View file

@ -0,0 +1,153 @@
#include "random_lcg.h"
#include <cmath>
// 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;
}

223
src/random_lcg.F90
View file

@ -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

56
src/random_lcg.h Normal file
View file

@ -0,0 +1,56 @@
#ifndef RANDOM_LCG_H
#define RANDOM_LCG_H
#include <cstdint>
//==============================================================================
// 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