Skip to content

Commit

Permalink
Merge branch 'mrnorman/cam-se/perf-sensitivity-merge' (PR #56)
Browse files Browse the repository at this point in the history
This branch contains two commits that add (1)a new method of seeding
Pseudo Random Numbers (PRN) used to perturb temperature and (2) a new
PRN Generator (PRNG) adapted by Pat Worley. The commits are
self-explanatory. These changes alter nothing in ACME by default. The
user must set new_random = .true. in the user_nl_cam in order to
activate the new code. I have specified file-by-file the changes
included in that file.

models/atm/cam/bld/namelist_files/namelist_definition.xml was changed
to add the new namelist parameters so that errors aren't thrown when
they're placed in the user_nl_cam.

models/atm/cam/src/control/cam_control_mod.F90 was changed to provide
a home for the new namelist variables.

models/atm/cam/src/control/runtime_opts.F90 was changed to document,
read in from namelist, and broadcast the new namelist variables to all
MPI tasks

models/atm/cam/src/dynamics/se/inidat.F90 was changed to allow the new
seeding process (if the new_random flag is true) and to use the new
PRNG (again if new_random == .true.). The seed_custom is then merged
with the existing element column GlobalID seed using a bitwise XOR,
and if seed_clock is .true. then the system_clock is also merged into
the seed with a bitwise XOR. This enables distinct seeds from one
ensemble run to the next for truly random numbers. Default behavior
remains identical to before and is only activated if the new flags
(which default to .false.) are set to .true.

models/atm/cam/src/dynamics/se/random_xgc.F90 is added to house the
new PRNG, which provides consistency across compilers and greater
fidelity since we don't know what's going on in most compilers'
PRNGs for the Fortran intrinsic.

[BFB]
  • Loading branch information
singhbalwinder committed Dec 18, 2014
2 parents c5a2860 + fa43c00 commit db4e72c
Show file tree
Hide file tree
Showing 5 changed files with 242 additions and 8 deletions.
21 changes: 21 additions & 0 deletions models/atm/cam/bld/namelist_files/namelist_definition.xml
Original file line number Diff line number Diff line change
Expand Up @@ -1318,6 +1318,27 @@ amount. Only applied for initial simulations.
Default: 0.0
</entry>

<entry id="new_random" type="logical" category="initial_conditions"
group="cam_inparm" valid_values="" >
Use a new RNG in dynamics/se/random_xgc.F90 that is more robust and repeatable than the compiler-dependent
fortran intrinsic random_number.
Default: FALSE
</entry>

<entry id="seed_custom" type="integer" category="initial_conditions"
group="cam_inparm" valid_values="" >
Specify a custom seed value for RNG when pertlim /= 0.0. If seed_custom == 0,
Then the default seed is used, which is the same from run to run, based only
on a unique column index
Default: 0
</entry>

<entry id="seed_clock" type="logical" category="initial_conditions"
group="cam_inparm" valid_values="" >
bitwise XOR the system_clock with the existing seed value to ensure it is unique.
Default: FALSE
</entry>

<entry id="readtrace" type="logical" category="initial_conditions"
group="cam_inparm" valid_values="" >
If TRUE, try to initialize data for all consituents by reading from the
Expand Down
5 changes: 4 additions & 1 deletion models/atm/cam/src/control/cam_control_mod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,10 @@ module cam_control_mod
!------------------------------------------------------------

! from perturb.h
real(r8) :: pertlim = 0.0_r8
real(r8) :: pertlim = 0.0_r8
logical :: new_random = .false.
integer :: seed_custom = 0
logical :: seed_clock = .false.

! from comadj.h
integer :: nlvdry = 3
Expand Down
17 changes: 16 additions & 1 deletion models/atm/cam/src/control/runtime_opts.F90
Original file line number Diff line number Diff line change
Expand Up @@ -194,6 +194,15 @@ module runtime_opts
! pertlim = n.n Max size of perturbation to apply to initial
! temperature field.
!
! new_random logical: if .true., use RNG in dynamics/se/random_xgc.F90
! instead of the fortran intrinsic.
!
! seed_custom integer: if > 0, use new seeding mechanism that uses a
! custom seed rather than a custom limit. Default 0
!
! seed_clock logical: if .true., XOR the system_clock with the seed,
! wheter it includes a custom seed or not. Default .false.
!
! phys_alltoall Dynamics/physics transpose option. See phys_grid module.
!
integer :: phys_alltoall
Expand Down Expand Up @@ -439,6 +448,9 @@ subroutine read_namelist(single_column_in, scmlon_in, scmlat_in, nlfilename_in )
dtime, &
nlvdry, &
pertlim ,&
new_random ,&
seed_custom ,&
seed_clock ,&
readtrace, rayk0, raykrange, raytau0, &
tracers_flag, &
inithist, indirect, &
Expand Down Expand Up @@ -900,7 +912,10 @@ subroutine distnl
call mpibcast (use_64bit_nc,1,mpilog,0,mpicom)
call mpibcast (print_step_cost,1,mpilog,0,mpicom)
call mpibcast (inithist_all ,1,mpilog,0,mpicom)
call mpibcast (pertlim ,1, mpir8, 0, mpicom )
call mpibcast (pertlim ,1, mpir8 , 0, mpicom )
call mpibcast (new_random ,1, mpilog, 0, mpicom )
call mpibcast (seed_custom ,1, mpiint, 0, mpicom )
call mpibcast (seed_clock ,1, mpilog, 0, mpicom )

call mpibcast (caseid ,len(caseid) ,mpichar,0,mpicom)
call mpibcast (avgflag_pertape, ptapes, mpichar,0,mpicom)
Expand Down
31 changes: 25 additions & 6 deletions models/atm/cam/src/dynamics/se/inidat.F90
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,8 @@ module inidat
use element_mod, only : element_t
use shr_kind_mod, only: r8 => shr_kind_r8
use spmd_utils, only: iam, masterproc
use cam_control_mod, only : ideal_phys, aqua_planet, pertlim
use cam_control_mod, only : ideal_phys, aqua_planet, pertlim, seed_custom, seed_clock, new_random
use random_xgc, only: init_ranx, ranx
implicit none
private
public read_inidat
Expand Down Expand Up @@ -71,6 +72,7 @@ subroutine read_inidat( ncid_ini, ncid_topo, dyn_in)
integer :: rndm_seed_sz
integer, allocatable :: rndm_seed(:)
real(r8) :: pertval
integer :: sysclk
integer :: i
real(r8), parameter :: D0_0 = 0.0_r8
real(r8), parameter :: D0_5 = 0.5_r8
Expand All @@ -86,7 +88,7 @@ subroutine read_inidat( ncid_ini, ncid_topo, dyn_in)

call get_dyn_decomp(elem, nlev, pio_double, iodesc)

lsize = pio_get_local_array_size(iodesc)
lsize = pio_get_local_array_size(iodesc)

tlncols = lsize/nlev

Expand Down Expand Up @@ -144,18 +146,35 @@ subroutine read_inidat( ncid_ini, ncid_topo, dyn_in)
'by +/- ', pertlim, ' to initial temperature field'
end if

call random_seed(size=rndm_seed_sz)
if (new_random) then
rndm_seed_sz = 1
else
call random_seed(size=rndm_seed_sz)
endif
allocate(rndm_seed(rndm_seed_sz))

do ie=1,nelemd
! seed random number generator based on element ID
! (possibly include a flag to allow clock-based random seeding)
rndm_seed = elem(ie)%GlobalId
call random_seed(put=rndm_seed)
rndm_seed(:) = elem(ie)%GlobalId
if (seed_custom > 0) rndm_seed(:) = ieor( rndm_seed(1) , int(seed_custom,kind(rndm_seed(1))) )
if (seed_clock) then
call system_clock(sysclk)
rndm_seed(:) = ieor( sysclk , int(rndm_seed(1),kind(sysclk)) )
endif
if (new_random) then
call init_ranx(rndm_seed(1))
else
call random_seed(put=rndm_seed)
endif
do i=1,np
do j=1,np
do k=1,nlev
call random_number(pertval)
if (new_random) then
pertval = ranx()
else
call random_number(pertval)
endif
pertval = D2_0*pertlim*(D0_5 - pertval)
elem(ie)%state%T(i,j,k,1) = elem(ie)%state%T(i,j,k,1)*(D1_0 + pertval)
end do
Expand Down
176 changes: 176 additions & 0 deletions models/atm/cam/src/dynamics/se/random_xgc.F90
Original file line number Diff line number Diff line change
@@ -0,0 +1,176 @@
! random number generator from XGC1 code
! imported November 19, 2014 by P. Worley (worleyph@ornl.gov)

MODULE Ecuyer_random
! L'Ecuyer's 1996 random number generator.
! Fortran version by Alan.Miller @ vic.cmis.csiro.au
! N.B. This version is compatible with Lahey's ELF90
! http://www.ozemail.com.au/~milleraj
! Latest revision - 30 March 1999

IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60)

! These are unsigned integers in the C version
INTEGER, SAVE :: s1 = 1234, s2 = -4567, s3 = 7890

CONTAINS

SUBROUTINE init_seeds(i1, i2, i3)
IMPLICIT NONE

INTEGER, INTENT(IN) :: i1, i2, i3

s1 = i1
s2 = i2
s3 = i3
IF (IAND(s1,-2) == 0) s1 = i1 - 1023
IF (IAND(s2,-8) == 0) s2 = i2 - 1023
IF (IAND(s3,-16) == 0) s3 = i3 - 1023

RETURN
END SUBROUTINE init_seeds

FUNCTION taus88() RESULT(random_numb)
! Generates a random number between 0 and 1. Translated from C function in:
! Reference:
! L'Ecuyer, P. (1996) `Maximally equidistributed combined Tausworthe
! generators', Math. of Comput., 65, 203-213.

! The cycle length is claimed to be about 2^(88) or about 3E+26.
! Actually - (2^31 - 1).(2^29 - 1).(2^28 - 1).

IMPLICIT NONE
REAL (dp) :: random_numb

INTEGER :: b

! N.B. ISHFT(i,j) is a bitwise (non-circular) shift operation;
! to the left if j > 0, otherwise to the right.

b = ISHFT( IEOR( ISHFT(s1,13), s1), -19)
s1 = IEOR( ISHFT( IAND(s1,-2), 12), b)
b = ISHFT( IEOR( ISHFT(s2,2), s2), -25)
s2 = IEOR( ISHFT( IAND(s2,-8), 4), b)
b = ISHFT( IEOR( ISHFT(s3,3), s3), -11)
s3 = IEOR( ISHFT( IAND(s3,-16), 17), b)
random_numb = IEOR( IEOR(s1,s2), s3) * 2.3283064365E-10_dp + 0.5_dp

RETURN
END FUNCTION taus88

SUBROUTINE init_seeds_ext(sv)
IMPLICIT NONE

INTEGER, INTENT(INOUT) :: sv(3)

IF (IAND(sv(1),-2) == 0) sv(1) = sv(1) - 1023
IF (IAND(sv(2),-8) == 0) sv(2) = sv(2) - 1023
IF (IAND(sv(3),-16) == 0) sv(3) = sv(3) - 1023

RETURN
END SUBROUTINE init_seeds_ext

FUNCTION taus88_ext(sv) RESULT(random_numb)
! Generates a random number between 0 and 1. Translated from C function in:
! Reference:
! L'Ecuyer, P. (1996) `Maximally equidistributed combined Tausworthe
! generators', Math. of Comput., 65, 203-213.

! The cycle length is claimed to be about 2^(88) or about 3E+26.
! Actually - (2^31 - 1).(2^29 - 1).(2^28 - 1).

IMPLICIT NONE
INTEGER, INTENT(INOUT) :: sv(3)
REAL (dp) :: random_numb

INTEGER :: b, i1, i2, i3

! N.B. ISHFT(i,j) is a bitwise (non-circular) shift operation;
! to the left if j > 0, otherwise to the right.

i1 = sv(1)
i2 = sv(2)
i3 = sv(3)
b = ISHFT( IEOR( ISHFT(i1,13), i1), -19)
i1 = IEOR( ISHFT( IAND(i1,-2), 12), b)
b = ISHFT( IEOR( ISHFT(i2,2), i2), -25)
i2 = IEOR( ISHFT( IAND(i2,-8), 4), b)
b = ISHFT( IEOR( ISHFT(i3,3), i3), -11)
i3 = IEOR( ISHFT( IAND(i3,-16), 17), b)
random_numb = IEOR( IEOR(i1,i2), i3) * 2.3283064365E-10_dp + 0.5_dp
sv(1) = i1
sv(2) = i2
sv(3) = i3

RETURN
END FUNCTION taus88_ext

END MODULE Ecuyer_random

!! random number generator module
module random_xgc
use Ecuyer_random, only: init_seeds, init_seeds_ext, taus88, taus88_ext
use shr_kind_mod, only: r8 => shr_kind_r8
implicit none
type seeds_type
integer, pointer :: s(:)
end type seeds_type
type(seeds_type), allocatable :: sv(:)

contains
!random number - seed initialize
! if ranx is to be executed in a threaded region,
! then each thread must set its own seed, otherwise
! each thread generates the identical sequence
subroutine init_ranx(seed)
implicit none

integer, intent(in) :: seed
#ifdef _OPENMP
integer :: omp_get_num_threads
integer :: omp_get_thread_num
integer :: num_threads, thread_id
integer :: i

thread_id = omp_get_thread_num()

!$OMP CRITICAL
num_threads = omp_get_num_threads()
if (.not. allocated(sv)) then
allocate( sv(0:num_threads-1) )
do i=0,num_threads-1
allocate(sv(i)%s(3))
sv(i)%s(1) = 1234*seed
sv(i)%s(2) = 2345*seed + 6789
sv(i)%s(3) = 4321*seed + 10
call init_seeds_ext(sv(i)%s)
end do
endif
!$OMP END CRITICAL

sv(thread_id)%s(1) = 1234*seed
sv(thread_id)%s(2) = 2345*seed + 6789
sv(thread_id)%s(3) = 4321*seed + 10
call init_seeds_ext(sv(thread_id)%s)

#else
call init_seeds(1234*seed,2345*seed+6789,4321*seed+10)
#endif
end subroutine init_ranx

! random number generator
function ranx()
implicit none
real(r8) :: ranx
#ifdef _OPENMP
integer :: thread_id
integer :: omp_get_thread_num
thread_id = omp_get_thread_num()
ranx=taus88_ext( sv(thread_id)%s )
#else
ranx=taus88()
#endif
end function ranx

end module random_xgc

0 comments on commit db4e72c

Please sign in to comment.