@brief Subroutine that samples a Gamma distribution representing the runaways' (marginal) energy distribution function. @details This subroutine uses the Metropolis-Hastings method for sampling the Gamma distribution representing the runaways' (marginal) energy distribution function. Unlike the typical Metropolis-Hasting method, after setting the boundaries of the region we want to sample, we perform a sampling in a larger region that contains the original sampling area plus a buffer region. After finishing the first sampling, we only keep the particles in the original sampling region, the particles in the p_buffer are sampled again until all of them lie within the original sampling region. This method ensures that the boundaries are well sampled.
@param[in] params Core KORC simulation parameters. @param[in,out] g Relativistic gamma factor @f$\gamma@f$ of the particles in a given species in the simulation. These are so that, they follow a Gamma distribution in energy. The parameters of the Gamma distributions are given by the user. @param[out] go Mean value of @f$\gamma@f$ of the particles in a given species. used to calculate the minimum required time step to resolve in detail the full-orbit dynamics of the particles. @param p Sampled normalized momentum @f$p' = p/m_ec@f$ of particles in a given particle species. @param p_buffer Size along the momentum axis of the buffer used in the Metropolis-Hastings method. @param p_test The test value of the normalized momentum used in the Metropolis-Hastings method. @param ratio Ratio of probabilities used to determine when a move in during the sampling is kept as part of the sampled chain. @param rand_unif A deviate of a uniform random distribution in the interval @f$[0,1]@f$. @param p_samples Temporary array to keep the sampled normalized momentum. @param deta Step size along the pitch-angle direction of the random walk used in the Metropolis-Hastings sampling. @param dp Step size along the momentum direction of the random walk used in the Metropolis-Hastings sampling. @param ii Iterator. @param ppp Number of particles per MPI processes. @param nsamples Number of total samples in the initial condition of a given simulation. @param mpierr MPI error status.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(KORC_PARAMS), | intent(in) | :: | params | |||
| real(kind=rp), | intent(inout), | DIMENSION(:), ALLOCATABLE | :: | g | ||
| real(kind=rp), | intent(out) | :: | go |
SUBROUTINE sample_gamma_distribution(params,g,go)
TYPE(KORC_PARAMS), INTENT(IN) :: params
REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: g
REAL(rp), INTENT(OUT) :: go
REAL(rp), DIMENSION(:), ALLOCATABLE :: p
REAL(rp) :: p_buffer
REAL(rp) :: p_test
REAL(rp) :: ratio
REAL(rp) :: rand_unif
REAL(rp), DIMENSION(:), ALLOCATABLE :: p_samples
REAL(rp) :: deta
REAL(rp) :: dp
INTEGER :: ii
INTEGER :: ppp
INTEGER :: nsamples
INTEGER :: mpierr
ppp = SIZE(g)
nsamples = ppp*params%mpi_params%nmpi
ALLOCATE(p(ppp))
dp = 1.0_rp
if (params%mpi_params%rank.EQ.0_idef) then
ALLOCATE(p_samples(nsamples))! Number of samples to distribute among all MPI processes
call RANDOM_SEED()
call RANDOM_NUMBER(rand_unif)
p_buffer = gamma_pdf_params%min_p + (gamma_pdf_params%max_p - gamma_pdf_params%min_p)*rand_unif
ii=2_idef
do while (ii .LE. 1000000_idef)
p_test = p_buffer + random_norm(0.0_rp,dp)
do while ((p_test.LT.gamma_pdf_params%min_p).OR.(p_test.GT.gamma_pdf_params%max_p))
p_test = p_buffer + random_norm(0.0_rp,dp)
end do
ratio = fRE(p_test)/fRE(p_buffer)
if (ratio .GE. 1.0_rp) then
p_buffer = p_test
ii = ii + 1_idef
else
call RANDOM_NUMBER(rand_unif)
if (rand_unif .LT. ratio) then
p_buffer = p_test
ii = ii + 1_idef
end if
end if
end do
call RANDOM_SEED()
call RANDOM_NUMBER(rand_unif)
p_samples(1) = p_buffer
ii=2_idef
do while (ii .LE. nsamples)
p_test = p_samples(ii-1) + random_norm(0.0_rp,dp)
do while ((p_test.LT.gamma_pdf_params%min_p).OR.(p_test.GT.gamma_pdf_params%max_p))
p_test = p_samples(ii-1) + random_norm(0.0_rp,dp)
end do
ratio = fRE(p_test)/fRE(p_samples(ii-1))
if (ratio .GE. 1.0_rp) then
p_samples(ii) = p_test
ii = ii + 1_idef
else
call RANDOM_NUMBER(rand_unif)
if (rand_unif .LT. ratio) then
p_samples(ii) = p_test
ii = ii + 1_idef
end if
end if
end do
go = SUM(SQRT(1.0_rp + p_samples**2))/nsamples
end if
CALL MPI_SCATTER(p_samples,ppp,MPI_REAL8,p,ppp,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
CALL MPI_BCAST(go,1,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
call MPI_BARRIER(MPI_COMM_WORLD,mpierr)
g = SQRT(1.0_rp + p**2)
DEALLOCATE(p)
if (params%mpi_params%rank.EQ.0_idef) then
DEALLOCATE(p_samples)
end if
END SUBROUTINE sample_gamma_distribution