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(inout) | :: | 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. |
||
| type(FIELDS), | intent(in) | :: | F | An instance of the KORC derived type FIELDS. |
subroutine sample_Hollmann_distribution_1Dtransport(params,spp,F)
!! @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(INOUT) :: 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.
TYPE(FIELDS), INTENT(IN) :: F
!! An instance of the KORC derived type FIELDS.
REAL(rp), DIMENSION(:), ALLOCATABLE :: R_samples,X_samples,Y_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 :: G_samples
!! Gamma of all samples
REAL(rp), DIMENSION(:), ALLOCATABLE :: eta_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) :: min_g,max_g
REAL(rp) :: min_R,max_R
REAL(rp) :: min_Z,max_Z
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) :: eta_buffer
!! Previous sample of pitch
REAL(rp) :: G_buffer
REAL(rp) :: R_test
!! Present sample of R location
REAL(rp) :: Z_test
!! Present sample of Z location
REAL(rp) :: eta_test
!! Present sample of pitch angle
REAL(rp) :: G_test
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) :: PSIp_lim,PSIP0,PSIN,PSIN0,PSIN1,sigma
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,rr
!! Sample iterator.
INTEGER :: mpierr
!! mpi error indicator
REAL(rp) :: dgmin,dgmax,deta
LOGICAL :: accepted
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/)
REAL(rp) :: EPHI,fRE_out,nAr0,nAr1,nAr2,nAr3,nD,nD1,ne,Te,Zeff,nRE
nsamples = spp%ppp*params%mpi_params%nmpi
psi_max_buff = spp%psi_max*1.25_rp
theta_rad=C_PI*spp%theta_gauss/180.0_rp
params%GC_coords=.TRUE.
PSIp_lim=F%PSIp_lim
PSIp0=F%PSIP_min
min_R=minval(F%X%R)
max_R=maxval(F%X%R)
min_Z=minval(F%X%Z)
max_Z=maxval(F%X%Z)
sigma=spp%sigmaR*params%cpp%length
!write(output_unit_write,*) min_R,max_R
!write(output_unit_write,*) min_Z,max_Z
deta = (h_params%max_pitch_angle - h_params%min_pitch_angle)/100.0_rp
dgmin = (h_params%min_sampling_g - h_params%min_g)/100.0_rp
dgmax = (h_params%max_g - h_params%max_sampling_g)/100.0_rp
if (h_params%min_sampling_g.gt.h_params%min_g) then
min_g=h_params%min_sampling_g
else
min_g=h_params%min_g
endif
max_g=h_params%max_sampling_g
! buffer at minimum gamma boundary
do ii=1_idef,INT(minmax_buffer_size,idef)
minmax = h_params%min_sampling_g - REAL(ii,rp)*dgmin
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(ii,rp)*dgmax
if (minmax.LT.h_params%max_g) then
max_g = minmax
end if
end do
!write(6,*) 'h_params%min_sampling_g',h_params%min_sampling_g,'h_params%min_g',h_params%min_g,'min_g',min_g
! buffer at minimum pitch angle boundary
if (h_params%min_pitch_angle.GE.korc_zero) then
do ii=1_idef,INT(minmax_buffer_size,idef)
minmax = h_params%min_pitch_angle - REAL(ii,rp)*deta
if (minmax.GT.0.0_rp) then
min_pitch_angle = minmax
end if
end do
else
min_pitch_angle = h_params%min_pitch_angle
end if
! buffer at maximum pitch angle boundary
do ii=1_idef,INT(minmax_buffer_size,idef)
minmax = h_params%max_pitch_angle + REAL(ii,rp)*deta
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))
ALLOCATE(X_samples(nsamples))
ALLOCATE(Y_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(eta_samples(nsamples))
! Number of samples to distribute among all MPI processes
ALLOCATE(G_samples(nsamples))
! Number of samples to distribute among all MPI processes
! Transient !
R_buffer = spp%Ro
Z_buffer = spp%Zo
if (.not.params%SameRandSeed) then
call init_random_seed()
else
call random_seed(put=seed)
end if
! initialize 2D gaussian argument and distribution function, or
! copy from previous sample
fRE_out=0._rp
do while (fRE_out.eq.0._rp)
eta_buffer = min_pitch_angle + (max_pitch_angle &
- min_pitch_angle)*get_random_U()
G_buffer = min_g + (max_g - min_g)*get_random_U()
!write(6,*) 'R_buffer',R_buffer*params%cpp%length
!write(6,*) 'Z_buffer',Z_buffer*params%cpp%length
!write(6,*) 'eta_buffer',eta_buffer
!write(6,*) 'G_buffer',G_buffer
call interp_nRE(params,R_buffer,0._rp,Z_buffer,psi0,EPHI,ne,Te,nRE, &
nAr0,nAr1,nAr2,nAr3,nD,nD1,G_buffer,fRE_out, &
rho1D=h_params%rho_axis(h_params%rho_ind))
!write(6,*) 'after first interp_nRE'
!call get_fields(params,spp%vars,F)
!psi0=spp%vars%PSI_P(1)
PSIN0=(psi0-PSIP0)/(PSIp_lim-PSIP0)
f0=nRE*fRE_out* &
fRE_H_pitch(params,eta_buffer,G_buffer,EPHI,ne,Te, &
nAr0,nAr1,nAr2,nAr3,nD,nD1)
!write(6,*) 'nRE',nRE*params%cpp%density
!write(6,*) 'fRE_out',fRE_out
end do
accepted=.false.
ii=1_idef
rr=1_idef
do while (ii .LE. 1000_idef)
if (modulo(ii,100).eq.0) then
write(output_unit_write,'("Burn: ",I10)') ii
write(6,'("Burn: ",I10)') ii
end if
R_test = R_buffer + get_random_N()*spp%dR
Z_test = Z_buffer + get_random_N()*spp%dZ
eta_test = eta_buffer + get_random_N()*spp%dth
G_test = G_buffer + get_random_N()*spp%dgam
! Test that pitch angle and momentum are within chosen boundary
do while ((eta_test .GT. max_pitch_angle).OR. &
(eta_test .LT. min_pitch_angle))
eta_test = eta_buffer + get_random_N()*spp%dth
end do
do while ((R_test.GT.max_R).OR.(R_test .LT. min_R))
R_test = R_buffer + get_random_N()*spp%dR
end do
do while ((Z_test.GT.max_Z).OR.(Z_test .LT. min_Z))
Z_test = Z_buffer + get_random_N()*spp%dZ
end do
do while ((G_test.LT.min_g).OR.(G_test.GT.max_g))
G_test = G_buffer + get_random_N()*spp%dgam
end do
if (accepted) then
!psi0=psi1
PSIN0=PSIN1
f0=f1
end if
! psi1=PSI_ROT_exp(R_test,spp%Ro,spp%sigmaR,Z_test,spp%Zo, &
! spp%sigmaZ,theta_rad)
!write(6,*) 'R_test',R_test*params%cpp%length
!write(6,*) 'Z_test',Z_test*params%cpp%length
!write(6,*) 'eta_test',eta_test
!write(6,*) 'G_test',G_test
call interp_nRE(params,R_test,0._rp,Z_test,psi1,EPHI,ne,Te,nRE, &
nAr0,nAr1,nAr2,nAr3,nD,nD1,G_test,fRE_out, &
rho1D=h_params%rho_axis(h_params%rho_ind))
PSIN1=(psi1-PSIP0)/(PSIp_lim-PSIP0)
!write(6,*) 'R_test',R_test*params%cpp%length
!write(6,*) 'Z_test',Z_test*params%cpp%length
!write(output_unit_write,*) 'ER',spp%vars%E(1,1)*params%cpp%Eo
!write(output_unit_write,*) 'EPHI',spp%vars%E(1,2)*params%cpp%Eo
!write(output_unit_write,*) 'EZ',spp%vars%E(1,3)*params%cpp%Eo
!write(output_unit_write,*) 'PSI',psi1
!write(output_unit_write,*) 'PSIN',PSIN
f1=nRE*fRE_out* &
fRE_H_pitch(params,eta_test,G_test,EPHI,ne,Te, &
nAr0,nAr1,nAr2,nAr3,nD,nD1)
! f1=fRE_H(eta_test,G_test)
!write(6,*) 'nRE',nRE*params%cpp%density
!write(6,*) 'fRE_out',fRE_out
! write(output_unit_write,'("psi0: ",E17.10)') psi0
! write(output_unit_write,'("psi1: ",E17.10)') psi1
!write(6,'("f0: ",E17.10)') f0
!write(6,'("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_exp(PSIN,spp%psi_max)* &
! R_test*f1*sin(deg2rad(eta_test))/ &
! (R_buffer*f0*sin(deg2rad(eta_buffer)))
ratio = indicator_exp(PSIN1,spp%psi_max)* &
R_test*f1*sin(deg2rad(eta_test))/ &
(R_buffer*f0*sin(deg2rad(eta_buffer)))
! ratio = f1*sin(deg2rad(eta_test))/(f0*sin(deg2rad(eta_buffer)))
if (rr.gt.100) then
write(6,*) 'f0,f1,R,Z,gam,eta',f0,f1,R_test*params%cpp%length,Z_test*params%cpp%length,G_test,eta_test
if (rr.gt.500) stop
endif
accepted=.false.
if (ratio .GE. 1.0_rp) then
accepted=.true.
R_buffer = R_test
Z_buffer = Z_test
eta_buffer = eta_test
G_buffer = G_test
ii = ii + 1_idef
rr=1_idef
else
! call RANDOM_NUMBER(rand_unif)
! if (rand_unif .LT. ratio) then
!if (get_random_mkl_U() .LT. ratio) then
if (get_random_U() .LT. ratio) then
accepted=.true.
R_buffer = R_test
Z_buffer = Z_test
eta_buffer = eta_test
G_buffer = G_test
ii = ii + 1_idef
rr=1_idef
else
rr=rr+1
end if
end if
end do
! Transient !
ii=1_idef
rr=1_idef
do while (ii .LE. nsamples)
! write(output_unit_write,'("sample:",I15)') ii
if (modulo(ii,nsamples/100).eq.0) then
write(output_unit_write,'("Sample: ",I10)') ii
write(6,'("Sample: ",I10)') ii
end if
R_test = R_buffer + get_random_N()*spp%dR
Z_test = Z_buffer + get_random_N()*spp%dZ
eta_test = eta_buffer + get_random_N()*spp%dth
G_test = G_buffer + get_random_N()*spp%dgam
! Test that pitch angle and momentum are within chosen boundary
do while ((eta_test .GT. max_pitch_angle).OR. &
(eta_test .LT. min_pitch_angle))
eta_test = eta_buffer + get_random_N()*spp%dth
end do
do while ((G_test.LT.min_g).OR.(G_test.GT.max_g))
G_test = G_buffer + get_random_N()*spp%dgam
end do
do while ((R_test.GT.max_R).OR.(R_test.LT. min_R))
R_test = R_buffer + get_random_N()*spp%dR
end do
do while ((Z_test.GT.max_Z).OR.(Z_test.LT. min_Z))
Z_test = Z_buffer + get_random_N()*spp%dZ
end do
if (accepted) then
PSIN0=PSIN1
f0=f1
end if
! psi1=PSI_ROT_exp(R_test,spp%Ro,spp%sigmaR,Z_test,spp%Zo, &
! spp%sigmaZ,theta_rad)
call interp_nRE(params,R_test,0._rp,Z_test,psi1,EPHI,ne,Te,nRE, &
nAr0,nAr1,nAr2,nAr3,nD,nD1,g_test,fRE_out, &
rho1D=h_params%rho_axis(h_params%rho_ind))
PSIN1=(psi1-PSIP0)/(PSIp_lim-PSIP0)
! write(output_unit_write,'("R: ",E17.10)') R_test
! write(output_unit_write,'("R0: ",E17.10)') spp%Ro
! write(output_unit_write,'("sigma_R: ",E17.10)') spp%sigmaR
! write(output_unit_write,'("dR: ",E17.10)') spp%dR
! write(output_unit_write,'("N_dR: ",E17.10)') random_norm(0.0_rp,spp%dR)
! write(output_unit_write,'("Z: ",E17.10)') Z_test
! write(output_unit_write,'("Z0: ",E17.10)') spp%Zo
! write(output_unit_write,'("sigma_Z: ",E17.10)') spp%sigmaZ
! write(output_unit_write,'("dZ: ",E17.10)') spp%dZ
! write(output_unit_write,'("N_dR: ",Z17.10)') random_norm(0.0_rp,spp%dZ)
f1=nRE*fRE_out* &
fRE_H_pitch(params,eta_test,G_test,EPHI,ne,Te, &
nAr0,nAr1,nAr2,nAr3,nD,nD1)
! f1=fRE_H(eta_test,G_test)
! ratio = indicator_exp(PSIN,psi_max_buff)* &
! R_test*f1*sin(deg2rad(eta_test))/ &
! (R_buffer*f0*sin(deg2rad(eta_buffer)))
ratio = indicator_exp(PSIN1,spp%psi_max)* &
R_test*f1*sin(deg2rad(eta_test))/ &
(R_buffer*f0*sin(deg2rad(eta_buffer)))
! ratio = f1*sin(deg2rad(eta_test))/(f0*sin(deg2rad(eta_buffer)))
if (rr.gt.100) then
write(6,*) 'f0,f1,R,Z,gam,eta',f0,f1,R_test*params%cpp%length,Z_test*params%cpp%length,G_test,eta_test
if (rr.gt.500) stop
endif
accepted=.false.
if (ratio .GE. 1.0_rp) then
accepted=.true.
R_buffer = R_test
Z_buffer = Z_test
eta_buffer = eta_test
G_buffer = G_test
rr=1_idef
else
!call RANDOM_NUMBER(rand_unif)
!if (rand_unif .LT. ratio) then
!if (get_random_mkl_U() .LT. ratio) then
if (get_random_U() .LT. ratio) then
accepted=.true.
R_buffer = R_test
Z_buffer = Z_test
eta_buffer = eta_test
G_buffer = G_test
rr=1_idef
else
rr=rr+1
end if
end if
! write(output_unit_write,'("R: ",E17.10)') R_buffer
! write(output_unit_write,'("Z: ",E17.10)') Z_buffer
! 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_exp(PSIN1,spp%psi_max)).EQ.1).AND. &
(G_buffer.LE.h_params%max_sampling_g).AND. &
(G_buffer.GE.h_params%min_sampling_g).AND. &
(eta_buffer.LE.h_params%max_pitch_angle).AND. &
(eta_buffer.GE.h_params%min_pitch_angle).AND. &
ACCEPTED) THEN
R_samples(ii) = R_buffer
Z_samples(ii) = Z_buffer
eta_samples(ii) = eta_buffer
G_samples(ii) = G_buffer
! write(output_unit_write,*) 'RS',R_buffer
! Sample phi location uniformly
!call RANDOM_NUMBER(rand_unif)
!PHI_samples(ii) = 2.0_rp*C_PI*rand_unif
!PHI_samples(ii) = 2.0_rp*C_PI*get_random_mkl_U()
PHI_samples(ii) = 2.0_rp*C_PI*get_random_U()
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
X_samples=R_samples*cos(PHI_samples)
Y_samples=R_samples*sin(PHI_samples)
! write(output_unit_write,*) 'R_samples',R_samples
! write(output_unit_write,*) 'PHI_samples',PHI_samples
! write(output_unit_write,*) 'Z_samples',Z_samples
! write(output_unit_write,*) 'G_samples',G_samples
! write(output_unit_write,*) 'eta_samples',eta_samples
if (TRIM(h_params%current_direction) .EQ. 'PARALLEL') then
eta_samples = 180.0_rp - eta_samples
end if
end if
params%GC_coords=.FALSE.
call MPI_BARRIER(MPI_COMM_WORLD,mpierr)
CALL MPI_SCATTER(X_samples,spp%ppp,MPI_REAL8, &
spp%vars%X(:,1),spp%ppp,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
CALL MPI_SCATTER(Y_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(eta_samples,spp%ppp,MPI_REAL8, &
spp%vars%eta,spp%ppp,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
CALL MPI_SCATTER(G_samples,spp%ppp,MPI_REAL8, &
spp%vars%g,spp%ppp,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
call MPI_BARRIER(MPI_COMM_WORLD,mpierr)
! write(output_unit_write,*) params%mpi_params%rank,'varX',spp%vars%X(:,1)
! write(output_unit_write,*) params%mpi_params%rank,'varY',spp%vars%X(:,2)
! 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,*) params%mpi_params%rank,'varX',spp%vars%X(:,1)
! write(output_unit_write,*) params%mpi_params%rank,'varR',spp%vars%Y(:,1)
! write(output_unit_write,'("Y_R: ",E17.10)') spp%vars%Y(:,1)*params%cpp%length
! write(output_unit_write,'("Y_PHI: ",E17.10)') spp%vars%Y(:,1)*params%cpp%length
! write(output_unit_write,'("Y_Z: ",E17.10)') spp%vars%Y(:,3)*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(X_samples)
DEALLOCATE(Y_samples)
DEALLOCATE(Z_samples)
DEALLOCATE(PHI_samples)
DEALLOCATE(eta_samples)
DEALLOCATE(G_samples)
end if
end subroutine sample_Hollmann_distribution_1Dtransport