Subroutine that generates a 2D Gaussian distribution in an elliptic torus as the initial spatial condition of a given particle species in the simulation.
*rand_unif
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(KORC_PARAMS), | intent(in) | :: | params | Core KORC simulation parameters. |
||
| type(SPECIES), | intent(inout) | :: | spp | An instance of the derived type SPECIES containing all the parameters and simulation variables of the different species in the simulation. |
subroutine MH_gaussian_elliptic_torus(params,spp)
!! @note Subroutine that generates a 2D Gaussian distribution in an
!! elliptic torus as the initial spatial condition of a given particle
!! species in the simulation. @endnote
TYPE(KORC_PARAMS), INTENT(IN) :: params
!! Core KORC simulation parameters.
TYPE(SPECIES), INTENT(INOUT) :: spp
!! An instance of the derived type SPECIES containing all the parameters
!! and simulation variables of the different species in the simulation.
REAL(rp), DIMENSION(:), ALLOCATABLE :: R_samples
!!
REAL(rp), DIMENSION(:), ALLOCATABLE :: Z_samples
!!
REAL(rp), DIMENSION(:), ALLOCATABLE :: ZETA_samples
!!
REAL(rp) :: psi_max_buff
!!
REAL(rp) :: theta_rad
!!
REAL(rp) :: R_buffer
!!
REAL(rp) :: Z_buffer
!!
REAL(rp) :: R_test
!!
REAL(rp) :: Z_test
!!
REAL(rp) :: psi0
!!
REAL(rp) :: psi1
!!
REAL(rp) :: rand_unif
!!
REAL(rp) :: ratio
!!
INTEGER :: nsamples
!!
INTEGER :: ii
!! Particle iterator.
INTEGER :: mpierr
LOGICAL :: accepted
nsamples = spp%ppp*params%mpi_params%nmpi
psi_max_buff = spp%psi_max*1.1_rp
theta_rad=C_PI*spp%theta_gauss/180.0_rp
if (params%mpi_params%rank.EQ.0_idef) then
ALLOCATE(R_samples(nsamples))
! Number of samples to distribute among all MPI processes
ALLOCATE(Z_samples(nsamples))
! Number of samples to distribute among all MPI processes
ALLOCATE(ZETA_samples(nsamples))
! Number of samples to distribute among all MPI processes
! Transient !
R_buffer = spp%Ro
Z_buffer = spp%Zo
ii=2_idef
do while (ii .LE. 1000_idef)
R_test = R_buffer + random_norm(0.0_rp,spp%sigmaR)
Z_test = Z_buffer + random_norm(0.0_rp,spp%sigmaZ)
psi0=PSI_ROT(R_buffer,spp%Ro,spp%sigmaR,Z_buffer,spp%Zo, &
spp%sigmaZ,theta_rad)
psi1=PSI_ROT(R_test,spp%Ro,spp%sigmaR,Z_test,spp%Zo,spp%sigmaZ,theta_rad)
ratio = indicator(psi1,spp%psi_max)*R_test*EXP(-psi1)/(R_buffer*EXP(-psi0))
if (ratio .GE. 1.0_rp) then
R_buffer = R_test
Z_buffer = Z_test
ii = ii + 1_idef
else
call RANDOM_NUMBER(rand_unif)
if (rand_unif .LT. ratio) then
R_buffer = R_test
Z_buffer = Z_test
ii = ii + 1_idef
end if
end if
end do
! Transient !
ii=1_idef
do while (ii .LE. nsamples)
R_test = R_buffer + random_norm(0.0_rp,spp%sigmaR)
Z_test = Z_buffer + random_norm(0.0_rp,spp%sigmaZ)
psi0=PSI_ROT(R_buffer,spp%Ro,spp%sigmaR,Z_buffer,spp%Zo, &
spp%sigmaZ,theta_rad)
psi1=PSI_ROT(R_test,spp%Ro,spp%sigmaR,Z_test,spp%Zo,spp%sigmaZ,theta_rad)
ratio = indicator(psi1,psi_max_buff)*R_test*EXP(-psi1)/(R_buffer*EXP(-psi0))
accepted=.false.
if (ratio .GE. 1.0_rp) then
accepted=.true.
R_buffer = R_test
Z_buffer = Z_test
else
call RANDOM_NUMBER(rand_unif)
if (rand_unif .LT. ratio) then
accepted=.true.
R_buffer = R_test
Z_buffer = Z_test
end if
end if
IF (INT(indicator(psi1,spp%psi_max)).EQ.1.and.accepted) THEN
R_samples(ii) = R_buffer
Z_samples(ii) = Z_buffer
! call RANDOM_NUMBER(rand_unif)
ZETA_samples(ii) = 2.0_rp*C_PI!!*rand_unif
ii = ii + 1_idef
END IF
end do
end if
CALL MPI_SCATTER(R_samples*sin(ZETA_samples),spp%ppp,MPI_REAL8, &
spp%vars%X(:,1),spp%ppp,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
CALL MPI_SCATTER(R_samples*cos(ZETA_samples),spp%ppp,MPI_REAL8, &
spp%vars%X(:,2),spp%ppp,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
CALL MPI_SCATTER(Z_samples,spp%ppp,MPI_REAL8,spp%vars%X(:,3), &
spp%ppp,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
call MPI_BARRIER(MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
DEALLOCATE(R_samples)
DEALLOCATE(Z_samples)
DEALLOCATE(ZETA_samples)
end if
end subroutine MH_gaussian_elliptic_torus