From 4d5ca78ba64c3494bf69a05698f299ea0c9c1bca Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 5 Oct 2017 12:33:12 -0500 Subject: [PATCH] Improve how PRN seed is exposed to Python --- openmc/capi/settings.py | 12 +++++ src/api.F90 | 7 +-- src/initialize.F90 | 6 ++- src/input_xml.F90 | 5 +- src/random_lcg.F90 | 109 ++++++++++++++++++++-------------------- 5 files changed, 77 insertions(+), 62 deletions(-) diff --git a/openmc/capi/settings.py b/openmc/capi/settings.py index 11b0eacef5..eede6cf47b 100644 --- a/openmc/capi/settings.py +++ b/openmc/capi/settings.py @@ -10,6 +10,10 @@ _RUN_MODES = {1: 'fixed source', 4: 'particle restart', 5: 'volume'} +_dll.openmc_set_seed.argtypes = [c_int64] +_dll.openmc_set_seed.restype = c_int +_dll.openmc_set_seed.errcheck = _error_handler + class _Settings(object): # Attributes that are accessed through a descriptor @@ -37,5 +41,13 @@ class _Settings(object): else: raise ValueError('Invalid run mode: {}'.format(mode)) + @property + def seed(self): + return c_int64.in_dll(_dll, 'seed').value + + @seed.setter + def seed(self, seed): + _dll.openmc_set_seed(seed) + settings = _Settings() diff --git a/src/api.F90 b/src/api.F90 index 3b2fdacad8..4410f2fe45 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -17,7 +17,7 @@ module openmc_api use initialize, only: openmc_init use particle_header, only: Particle use plot, only: openmc_plot_geometry - use random_lcg, only: seed, initialize_prng + use random_lcg, only: seed, openmc_set_seed use settings use simulation_header use tally_header @@ -212,6 +212,8 @@ contains !=============================================================================== subroutine openmc_hard_reset() bind(C) + integer :: err + ! Reset all tallies and timers call openmc_reset() @@ -220,8 +222,7 @@ contains total_gen = 0 ! Reset the random number generator state - seed = 1_8 - call initialize_prng() + err = openmc_set_seed(1_8) end subroutine openmc_hard_reset !=============================================================================== diff --git a/src/initialize.F90 b/src/initialize.F90 index 1720f2adf9..a5ead1e835 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -21,7 +21,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: initialize_prng + use random_lcg, only: openmc_set_seed, seed use settings #ifdef _OPENMP use simulation_header, only: n_threads @@ -45,6 +45,8 @@ contains subroutine openmc_init(intracomm) bind(C) integer, intent(in), optional :: intracomm ! MPI intracommunicator + integer :: err + ! Copy the communicator to a new variable. This is done to avoid changing ! the signature of this subroutine. If MPI is being used but no communicator ! was passed, assume MPI_COMM_WORLD. @@ -83,7 +85,7 @@ contains ! Initialize random number generator -- if the user specifies a seed, it ! will be re-initialized later - call initialize_prng() + err = openmc_set_seed(seed) ! Read XML input files call read_input_xml() diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 545c7e3510..402e030c41 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -25,7 +25,7 @@ module input_xml use nuclide_header use output, only: write_message, title, header, print_plot use plot_header - use random_lcg, only: prn, seed, initialize_prng + use random_lcg, only: prn, openmc_set_seed use surface_header use set_header, only: SetChar use settings @@ -151,6 +151,7 @@ contains integer :: temp_int integer :: temp_int_array3(3) integer(C_INT32_T) :: i_start, i_end + integer(C_INT64_T) :: seed integer(C_INT) :: err integer, allocatable :: temp_int_array(:) integer :: n_tracks @@ -340,7 +341,7 @@ contains ! Copy random number seed if specified if (check_for_node(root, "seed")) then call get_node_value(root, "seed", seed) - call initialize_prng() + err = openmc_set_seed(seed) end if ! Number of bins for logarithmic grid diff --git a/src/random_lcg.F90 b/src/random_lcg.F90 index 287bbba6d9..685557f070 100644 --- a/src/random_lcg.F90 +++ b/src/random_lcg.F90 @@ -1,5 +1,7 @@ module random_lcg + use, intrinsic :: ISO_C_BINDING + use constants implicit none @@ -7,28 +9,29 @@ module random_lcg private save - ! Random number seed - integer(8), public :: seed = 1_8 + ! Starting seed + integer(C_INT64_T), public, bind(C) :: seed = 1_8 - integer(8) :: prn_seed0 ! original seed - integer(8) :: prn_seed(N_STREAMS) ! current seed - integer(8) :: prn_mult ! multiplication factor, g - integer(8) :: prn_add ! additive factor, c - integer :: prn_bits ! number of bits, M - integer(8) :: prn_mod ! 2^M - integer(8) :: prn_mask ! 2^M - 1 - integer(8) :: prn_stride ! stride between particles - real(8) :: prn_norm ! 2^(-M) - integer :: stream ! current RNG stream + ! 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) + ! Current PRNG state + integer(C_INT64_T) :: prn_seed(N_STREAMS) ! current seed + integer :: stream ! current RNG stream !$omp threadprivate(prn_seed, stream) public :: prn public :: future_prn - public :: initialize_prng public :: set_particle_seed public :: advance_prn_seed public :: prn_set_stream + public :: openmc_set_seed contains @@ -38,10 +41,10 @@ contains function prn() result(pseudo_rn) - real(8) :: pseudo_rn + real(C_DOUBLE) :: pseudo_rn - ! This algorithm uses bit-masking to find the next integer(8) value to be - ! used to calculate the random number + ! 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) @@ -59,40 +62,14 @@ contains function future_prn(n) result(pseudo_rn) - integer(8), intent(in) :: n ! number of prns to skip + integer(C_INT64_T), intent(in) :: n ! number of prns to skip - real(8) :: pseudo_rn + real(C_DOUBLE) :: pseudo_rn pseudo_rn = future_seed(n, prn_seed(stream)) * prn_norm end function future_prn -!=============================================================================== -! INITIALIZE_PRNG sets up the random number generator, determining the seed and -! values for g, c, and m. -!=============================================================================== - - subroutine initialize_prng() - - integer :: i - - prn_seed0 = seed -!$omp parallel - do i = 1, N_STREAMS - prn_seed(i) = prn_seed0 + i - 1 - end do - stream = STREAM_TRACKING -!$omp end parallel - prn_mult = 2806196910506780709_8 - prn_add = 1_8 - prn_bits = 63 - prn_mod = ibset(0_8, prn_bits) ! clever way of calculating 2**bits - prn_mask = prn_mod - 1_8 - prn_stride = 152917_8 - prn_norm = 2._8**(-prn_bits) - - end subroutine initialize_prng - !=============================================================================== ! SET_PARTICLE_SEED sets the seed to a unique value based on the ID of the ! particle @@ -100,12 +77,12 @@ contains subroutine set_particle_seed(id) - integer(8), intent(in) :: id + integer(C_INT64_T), intent(in) :: id integer :: i do i = 1, N_STREAMS - prn_seed(i) = future_seed(id*prn_stride, prn_seed0 + i - 1) + prn_seed(i) = future_seed(id*prn_stride, seed + i - 1) end do end subroutine set_particle_seed @@ -117,7 +94,7 @@ contains subroutine advance_prn_seed(n) - integer(8), intent(in) :: n ! number of seeds to skip + integer(C_INT64_T), intent(in) :: n ! number of seeds to skip prn_seed(stream) = future_seed(n, prn_seed(stream)) @@ -132,15 +109,15 @@ contains function future_seed(n, seed) result(new_seed) - integer(8), intent(in) :: n ! number of seeds to skip - integer(8), intent(in) :: seed ! original seed - integer(8) :: new_seed ! 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(8) :: nskip ! positive number of seeds to skip - integer(8) :: g ! original multiplicative constant - integer(8) :: c ! original additive constnat - integer(8) :: g_new ! new effective multiplicative constant - integer(8) :: c_new ! new effective additive constant + 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 @@ -198,4 +175,26 @@ contains 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 + end module random_lcg