mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 13:45:36 -04:00
Improve how PRN seed is exposed to Python
This commit is contained in:
Paul Romano
parent
12841f11af
commit
4d5ca78ba6
5 changed files with 77 additions and 62 deletions
|
|
@ -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()
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
!===============================================================================
|
||||
|
|
|
|||
|
|
@ -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()
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue