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

Merge pull request #939 from smharper/c_rng

Browse source 
Translate the RNG to C++
This commit is contained in:
Paul Romano 2017-12-15 10:13:13 +07:00 • committed by GitHub
commit 25ce61622d
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
6 changed files with 317 additions and 200 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 "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)
#===============================================================================

15
src/constants.F90
View file

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

4
src/initialize.F90
View file

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

222
src/random_lcg.F90
View file

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

145
src/random_lcg.cpp Normal file
View file

@ -0,0 +1,145 @@
#include "random_lcg.h"
#include <cmath>
// 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<uint64_t>(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<uint64_t>(id) * prn_stride, seed + i);
}
}
//==============================================================================
// ADVANCE_PRN_SEED
//==============================================================================
extern "C" void
advance_prn_seed(int64_t n)
{
prn_seed[stream] = future_seed(static_cast<uint64_t>(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;
}

85
src/random_lcg.h Normal file
View file

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