Subroutine that generates a 2D Gaussian distribution in an elliptic torus as the initial spatial condition of a given particle species in the simulation.
| 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. An instance of the KORC derived type FIELDS. |
subroutine Spong_3D(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.
!! An instance of the KORC derived type FIELDS.
REAL(rp), DIMENSION(:), ALLOCATABLE :: R_samples
!! Major radial location of all samples
REAL(rp), DIMENSION(:), ALLOCATABLE :: PHI_samples
!! Azimuithal angle of all samples
REAL(rp), DIMENSION(:), ALLOCATABLE :: Z_samples
!! Vertical location of all samples
REAL(rp), DIMENSION(:), ALLOCATABLE :: T_samples
!! Pitch angle of all samples
REAL(rp) :: psi_max_buff
!! Value of buffer above desired maximum argument of 2D Gaussian spatial
!! profile
REAL(rp) :: minmax
!! Temporary variable used for setting buffers
!! Minimum domain for momentum sampling including buffer
REAL(rp) :: max_pitch_angle
!! Maximum domain for pitch angle sampling including buffer
REAL(rp) :: min_pitch_angle
!! Minimum domain for pitch angle sampling including buffer
REAL(rp) :: theta_rad
!! Angle of rotation of 2D Gaussian spatial distribution in radians
REAL(rp) :: R_buffer
!! Previous sample of R location
REAL(rp) :: Z_buffer
!! Previous sample of Z location
REAL(rp) :: T_buffer
!! Previous sample of pitch angle
REAL(rp) :: R_test
!! Present sample of R location
REAL(rp) :: Z_test
!! Present sample of Z location
REAL(rp) :: T_test
!! Present sample of pitch angle
REAL(rp) :: psi0
!! Previous value of 2D Gaussian argument based on R_buffer, Z_buffer
REAL(rp) :: psi1
!! Present value of 2D Gaussian argument based on R_test, Z_test
REAL(rp) :: f0
!! Evaluation of Avalanche distribution with previous sample
REAL(rp) :: f1
!! Evaluation of Avalanche distribution with present sample
REAL(rp) :: rand_unif
!! Uniform random variable [0,1]
REAL(rp) :: ratio
!! MH selection criteria
INTEGER :: nsamples
!! Total number of samples to be distributed over all mpi processes
INTEGER :: ii
!! Sample iterator.
INTEGER :: mpierr
!! mpi error indicator
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
! buffer at minimum pitch angle boundary
if (spp%etao_lims(1).GE.korc_zero) then
do ii=1_idef,INT(minmax_buffer_size,idef)
minmax = spp%etao_lims(1) - REAL(ii,rp)* &
(spp%etao_lims(2)-spp%etao_lims(1))/100_rp
if (minmax.GT.0.0_rp) then
min_pitch_angle = minmax
end if
end do
else
min_pitch_angle = spp%etao_lims(1)
end if
! buffer at maximum pitch angle boundary
do ii=1_idef,INT(minmax_buffer_size,idef)
minmax = spp%etao_lims(2) + REAL(ii,rp)* &
(spp%etao_lims(2)-spp%etao_lims(1))/100_rp
if (minmax.LE.180.0_rp) then
max_pitch_angle = minmax
end if
end do
if (params%mpi_params%rank.EQ.0_idef) then
ALLOCATE(R_samples(nsamples))
! Number of samples to distribute among all MPI processes
ALLOCATE(PHI_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(T_samples(nsamples))
! Number of samples to distribute among all MPI processes
! Transient !
R_buffer = spp%Ro
Z_buffer = spp%Zo
call RANDOM_NUMBER(rand_unif)
T_buffer = min_pitch_angle + (max_pitch_angle &
- min_pitch_angle)*rand_unif
! write(output_unit_write,'("length norm: ",E17.10)') params%cpp%length
ii=1_idef
do while (ii .LE. 1000_idef)
! write(output_unit_write,'("burn:",I15)') ii
R_test = R_buffer + random_norm(0.0_rp,spp%dR)
Z_test = Z_buffer + random_norm(0.0_rp,spp%dZ)
T_test = T_buffer + random_norm(0.0_rp,spp%dth)
! Test that pitch angle and momentum are within chosen boundary
do while ((T_test .GT. spp%etao_lims(2)).OR. &
(T_test .LT. spp%etao_lims(1)))
T_test = T_buffer + random_norm(0.0_rp,spp%dth)
end do
! initialize 2D gaussian argument and distribution function, or
! copy from previous sample
if (ii==1) then
psi0=PSI_ROT(R_buffer,spp%Ro,spp%sigmaR,Z_buffer,spp%Zo, &
spp%sigmaZ,theta_rad)
f0=Spong_2D(spp%Ro,spp%Spong_b,spp%Spong_w,spp%Spong_dlam, &
R_buffer,Z_buffer,T_buffer)
else
psi0=psi1
f0=f1
end if
psi1=PSI_ROT(R_test,spp%Ro,spp%sigmaR,Z_test,spp%Zo, &
spp%sigmaZ,theta_rad)
f1=Spong_2D(spp%Ro,spp%Spong_b,spp%Spong_w,spp%Spong_dlam, &
R_test,Z_test,T_test)
! write(output_unit_write,'("psi0: ",E17.10)') psi0
! write(output_unit_write,'("psi1: ",E17.10)') psi1
! write(output_unit_write,'("f0: ",E17.10)') f0
! write(output_unit_write,'("f1: ",E17.10)') f1
! Calculate acceptance ratio for MH algorithm. fRE function
! incorporates p^2 factor of spherical coordinate Jacobian
! for velocity phase space, factors of sin(pitch angle) for velocity
! phase space and cylindrical coordinate Jacobian R for spatial
! phase space incorporated here.
ratio = indicator(psi1,spp%psi_max)*R_test*EXP(-psi1)*f1/ &
(R_buffer*EXP(-psi0)*f0)
if (ratio .GE. 1.0_rp) then
R_buffer = R_test
Z_buffer = Z_test
T_buffer = T_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
T_buffer = T_test
ii = ii + 1_idef
end if
end if
end do
! Transient !
ii=1_idef
do while (ii .LE. nsamples)
! write(output_unit_write,'("sample:",I15)') ii
if (modulo(ii,10000).eq.0) then
write(output_unit_write,'("Sample: ",I10)') ii
end if
R_test = R_buffer + random_norm(0.0_rp,spp%dR)
Z_test = Z_buffer + random_norm(0.0_rp,spp%dZ)
T_test = T_buffer + random_norm(0.0_rp,spp%dth)
! Selection boundary is set with buffer region
do while ((T_test .GT. max_pitch_angle).OR. &
(T_test .LT. min_pitch_angle))
if (T_test.lt.0) then
T_test=abs(T_test)
exit
end if
T_test = T_buffer + random_norm(0.0_rp,spp%dth)
end do
psi0=psi1
f0=f1
psi1=PSI_ROT(R_test,spp%Ro,spp%sigmaR,Z_test,spp%Zo, &
spp%sigmaZ,theta_rad)
f1=Spong_2D(spp%Ro,spp%Spong_b,spp%Spong_w,spp%Spong_dlam, &
R_test,Z_test,T_test)
ratio = indicator(psi1,psi_max_buff)*R_test*EXP(-psi1)*f1/ &
(R_buffer*EXP(-psi0)*f0)
if (ratio .GE. 1.0_rp) then
R_buffer = R_test
Z_buffer = Z_test
T_buffer = T_test
else
call RANDOM_NUMBER(rand_unif)
if (rand_unif .LT. ratio) then
R_buffer = R_test
Z_buffer = Z_test
T_buffer = T_test
end if
end if
! Only accept sample if it is within desired boundary, but
! add to MC above if within buffer. This helps make the boundary
! more defined.
IF ((INT(indicator(psi1,spp%psi_max)).EQ.1).AND. &
(T_buffer.LE.spp%etao_lims(2)).AND. &
(T_buffer.GE.spp%etao_lims(1))) THEN
R_samples(ii) = R_buffer
Z_samples(ii) = Z_buffer
T_samples(ii) = T_buffer
! Sample phi location uniformly
call RANDOM_NUMBER(rand_unif)
PHI_samples(ii) = 2.0_rp*C_PI*rand_unif
ii = ii + 1_idef
END IF
end do
! if (minval(R_samples(:)).lt.1._rp/params%cpp%length) stop 'error with sample'
! write(output_unit_write,'("R_sample: ",E17.10)') R_samples(:)*params%cpp%length
end if
CALL MPI_SCATTER(R_samples*cos(PHI_samples),spp%ppp,MPI_REAL8, &
spp%vars%X(:,1),spp%ppp,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
CALL MPI_SCATTER(R_samples*sin(PHI_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_SCATTER(T_samples,spp%ppp,MPI_REAL8, &
! spp%vars%eta,spp%ppp,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
call MPI_BARRIER(MPI_COMM_WORLD,mpierr)
! write(output_unit_write,'("X_X: ",E17.10)') spp%vars%X(:,1)*params%cpp%length
! gamma is kept for each particle, not the momentum
if (params%orbit_model(1:2).eq.'GC') call cart_to_cyl(spp%vars%X,spp%vars%Y)
! write(output_unit_write,'("Y_R: ",E17.10)') spp%vars%Y(:,1)*params%cpp%length
! if (minval(spp%vars%Y(:,1)).lt.1._rp/params%cpp%length) stop 'error with avalanche'
if (params%mpi_params%rank.EQ.0_idef) then
DEALLOCATE(R_samples)
DEALLOCATE(Z_samples)
DEALLOCATE(PHI_samples)
DEALLOCATE(T_samples)
end if
end subroutine Spong_3D