MCMC and MH algorithm perfomred on single MPI process to sample distribution function fRE_H
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(KORC_PARAMS), | intent(in) | :: | params | |||
| type(SPECIES), | intent(inout) | :: | spp |
SUBROUTINE sample_Hollmann_distribution(params,spp)
TYPE(KORC_PARAMS), INTENT(IN) :: params
TYPE(SPECIES), INTENT(INOUT) :: spp
! REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: g
! REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: eta
! REAL(rp), INTENT(OUT) :: go
! REAL(rp), INTENT(OUT) :: etao
REAL(rp), DIMENSION(:), ALLOCATABLE :: p
REAL(rp) :: g_buffer
REAL(rp) :: g_test
REAL(rp) :: eta_buffer
REAL(rp) :: eta_test
REAL(rp) :: ratio
REAL(rp) :: rand_unif
REAL(rp), DIMENSION(:), ALLOCATABLE :: g_samples
REAL(rp), DIMENSION(:), ALLOCATABLE :: eta_samples
REAL(rp), DIMENSION(:), ALLOCATABLE :: g_tmp
REAL(rp), DIMENSION(:), ALLOCATABLE :: eta_tmp
REAL(rp) :: minmax
REAL(rp) :: min_g
REAL(rp) :: max_g
REAL(rp) :: min_pitch_angle
REAL(rp) :: max_pitch_angle
REAL(rp) :: dg
REAL(rp) :: deta
LOGICAL :: lp
INTEGER :: index_i
INTEGER :: index_f
INTEGER :: num_accepted
INTEGER :: ii
INTEGER :: jj
INTEGER :: ppp
INTEGER :: nsamples
INTEGER :: mpierr
INTEGER,DIMENSION(33) :: seed=(/1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1/)
nsamples = spp%ppp*params%mpi_params%nmpi
index_i = MINLOC(ABS(h_params%g - h_params%min_sampling_g),1)
!index of minimum gamma range desired
index_f = MINLOC(ABS(h_params%g - h_params%max_sampling_g),1)
!index of maximum gamma range desired
deta = (h_params%max_pitch_angle - h_params%min_pitch_angle)/100.0_rp
dg = (h_params%max_sampling_g - h_params%min_sampling_g)/100.0_rp
! buffer at minimum gamma boundary
do jj=1_idef,INT(minmax_buffer_size,idef)
minmax = h_params%min_sampling_g - REAL(jj,rp)*dg
if (minmax.GT.h_params%min_g) then
min_g = minmax
end if
! buffer at maximum gamma boundary
minmax = h_params%max_sampling_g + REAL(jj,rp)*dg
if (minmax.LT.h_params%max_g) then
max_g = minmax
end if
end do
! buffer at minimum pitch angle boundary
if (h_params%min_pitch_angle.GE.korc_zero) then
do jj=1_idef,INT(minmax_buffer_size,idef)
minmax = h_params%min_pitch_angle - REAL(jj,rp)*deta
if (minmax.GT.0.0_rp) then
min_pitch_angle = minmax
end if
end do
else
min_pitch_angle = 0.0_rp
end if
! buffer at maximum pitch angle boundary
do jj=1_idef,INT(minmax_buffer_size,idef)
minmax = h_params%max_pitch_angle + REAL(jj,rp)*deta
if (minmax.LE.90.0_rp) then
max_pitch_angle = minmax
else
max_pitch_angle = h_params%max_pitch_angle
EXIT
end if
end do
if (params%mpi_params%rank.EQ.0_idef) then
!! MCMC and MH algorithm perfomred on single MPI process
!! to sample distribution function fRE_H
ALLOCATE(g_samples(nsamples))
! Number of samples to distribute among all MPI processes
ALLOCATE(eta_samples(nsamples))
! Number of samples to distribute among all MPI processes
ALLOCATE(g_tmp(nsamples))
! Number of samples to distribute among all MPI processes
ALLOCATE(eta_tmp(nsamples))
! Number of samples to distribute among all MPI processes
!Transient!
if (.not.params%SameRandSeed) then
call init_random_seed()
else
call random_seed(put=seed)
end if
call RANDOM_NUMBER(rand_unif)
! rand_unif=get_random_U()
eta_buffer = h_params%min_pitch_angle + (h_params%max_pitch_angle &
- h_params%min_pitch_angle)*rand_unif
call RANDOM_NUMBER(rand_unif)
! rand_unif=get_random_U()
g_buffer = h_params%min_sampling_g + (h_params%max_sampling_g - &
h_params%min_sampling_g)*rand_unif
ii=2_idef
do while (ii .LE. 1000_idef)
eta_test = eta_buffer + random_norm(0.0_rp,spp%dth)
! eta_test = eta_buffer + get_random_N()*spp%dth
do while ((ABS(eta_test) .GT. h_params%max_pitch_angle).OR. &
(ABS(eta_test) .LT. h_params%min_pitch_angle))
eta_test = eta_buffer + random_norm(0.0_rp,spp%dth)
! eta_test = eta_buffer + get_random_N()*spp%dth
end do
g_test = g_buffer + random_norm(0.0_rp,spp%dgam)
! g_test = g_buffer + get_random_N()*spp%dgam
do while ((g_test.LT.h_params%min_sampling_g).OR. &
(g_test.GT.h_params%max_sampling_g))
g_test = g_buffer + random_norm(0.0_rp,spp%dgam)
! g_test = g_buffer + get_random_N()*spp%dgam
end do
ratio = fRE_H(eta_test,g_test)*sin(deg2rad(eta_test))/ &
(fRE_H(eta_buffer,g_buffer)* &
sin(deg2rad(eta_buffer)))
!ratio = fRE_H(eta_test,g_test)/fRE_H(eta_buffer,g_buffer)
!ratio = fRE_HxPR(eta_test,g_test)/fRE_HxPR(eta_buffer,g_buffer)
if (ratio .GE. 1.0_rp) then
g_buffer = g_test
eta_buffer = eta_test
ii = ii + 1_idef
else
call RANDOM_NUMBER(rand_unif)
! rand_unif=get_random_U()
if (rand_unif .LT. ratio) then
g_buffer = g_test
eta_buffer = eta_test
ii = ii + 1_idef
end if
end if
end do
!Transient!
eta_tmp(1) = eta_buffer
g_tmp(1) = g_buffer
num_accepted = 0_idef
do while(num_accepted.LT.nsamples)
ii=2_idef
do while (ii .LE. nsamples)
if (modulo(ii,nsamples/10).eq.0) then
write(output_unit_write,'("Sample: ",I10)') ii
end if
! write(output_unit_write,'("iisample",I16)') ii
eta_test = eta_tmp(ii-1) + random_norm(0.0_rp,spp%dth)
!eta_test = eta_tmp(ii-1) + get_random_N()*spp%dth
! write(output_unit_write,'("max_pitch_angle: ",E17.10)') max_pitch_angle
! write(output_unit_write,'("min_pitch_angle: ",E17.10)') min_pitch_angle
do while ((ABS(eta_test) .GT. max_pitch_angle).OR. &
(ABS(eta_test) .LT. min_pitch_angle))
eta_test = eta_tmp(ii-1) + random_norm(0.0_rp,spp%dth)
! eta_test = eta_tmp(ii-1) + get_random_N()*spp%dth
! write(output_unit_write,'("eta_test: ",E17.10)') eta_test
end do
g_test = g_tmp(ii-1) + random_norm(0.0_rp,spp%dgam)
!g_test = g_tmp(ii-1) + get_random_N()*spp%dgam
! write(output_unit_write,'("max_g: ",E17.10)') max_g
! write(output_unit_write,'("min_g: ",E17.10)') min_g
do while ((g_test.LT.min_g).OR.(g_test.GT.max_g))
g_test = g_tmp(ii-1) + random_norm(0.0_rp,spp%dgam)
! g_test = g_tmp(ii-1) + get_random_N()*spp%dgam
! write(output_unit_write,'("g_test: ",E17.10)') g_test
end do
ratio = fRE_H(eta_test,g_test)*sin(deg2rad(eta_test))/ &
(fRE_H(eta_tmp(ii-1),g_tmp(ii-1))* &
sin(deg2rad(eta_tmp(ii-1))))
! ratio = fRE_H(eta_test,g_test)/fRE_H(eta_tmp(ii-1),g_tmp(ii-1))
!ratio = fRE_HxPR(eta_test,g_test)/fRE_HxPR(eta_tmp(ii-1),g_tmp(ii-1))
! write(output_unit_write,'("ratio: ",E17.10)') ratio
if (ratio .GE. 1.0_rp) then
g_tmp(ii) = g_test
eta_tmp(ii) = eta_test
ii = ii + 1_idef
else
call RANDOM_NUMBER(rand_unif)
! rand_unif=get_random_U()
if (rand_unif .LT. ratio) then
g_tmp(ii) = g_test
eta_tmp(ii) = eta_test
ii = ii + 1_idef
end if
end if
end do
eta_tmp = ABS(eta_tmp)
ii = 1_idef
do while ( (ii.LT.nsamples).AND.(num_accepted.LT.nsamples) )
! write(output_unit_write,'("iiaccept",I16)') ii
lp = (g_tmp(ii).LE.h_params%max_sampling_g).AND. &
(g_tmp(ii).GE.h_params%min_sampling_g).AND. &
(eta_tmp(ii).LE.h_params%max_pitch_angle).AND. &
(eta_tmp(ii).GE.h_params%min_pitch_angle)
if (lp) then
num_accepted = num_accepted + 1_idef
g_samples(num_accepted) = g_tmp(ii)
eta_samples(num_accepted) = eta_tmp(ii)
end if
ii = ii + 1_idef
end do
end do
! if (TRIM(h_params%current_direction) .EQ. 'ANTICLOCKWISE') then
! eta_samples = 180.0_rp - eta_samples
! end if
! go = SUM(g_samples)/nsamples
! etao = SUM(eta_samples)/nsamples
end if !MCMC computed on single MPI process
CALL MPI_SCATTER(g_samples,spp%ppp,MPI_REAL8,spp%vars%g,spp%ppp,MPI_REAL8, &
0,MPI_COMM_WORLD,mpierr)
CALL MPI_SCATTER(eta_samples,spp%ppp,MPI_REAL8,spp%vars%eta,spp%ppp,MPI_REAL8, &
0,MPI_COMM_WORLD,mpierr)
! CALL MPI_BCAST(go,1,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
! CALL MPI_BCAST(etao,1,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
call MPI_BARRIER(MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
DEALLOCATE(g_samples)
DEALLOCATE(eta_samples)
DEALLOCATE(g_tmp)
DEALLOCATE(eta_tmp)
end if
! write(output_unit_write,'("sampled eta: ",E17.10)') eta
END SUBROUTINE sample_Hollmann_distribution