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. |
||
| type(PROFILES), | intent(in) | :: | P | An instance of the KORC derived type PROFILES. |
||
| type(FIELDS), | intent(in) | :: | F | An instance of the KORC derived type FIELDS. |
subroutine Avalanche_4D(params,spp,P,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(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.
TYPE(PROFILES), INTENT(IN) :: P
!! An instance of the KORC derived type PROFILES.
TYPE(FIELDS), INTENT(IN) :: F
!! 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 :: P_samples
!! Magnitude of momentum of all samples
REAL(rp), DIMENSION(:), ALLOCATABLE :: T_samples
!! Pitch angle of all samples
REAL(rp), DIMENSION(:), ALLOCATABLE :: mom
!! Magnitude of momentum of samples on individual mpi process
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
REAL(rp) :: max_p
!! Maximum domain for momentum sampling including buffer
REAL(rp) :: min_p
!! 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) :: P_buffer
!! Previous sample of momentum
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) :: P_test
!! Present sample of momentum
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) :: fRE0
!! Evaluation of Avalanche distribution with previous sample
REAL(rp) :: fRE1
!! 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
ALLOCATE(mom(spp%ppp))
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 p boundary
do ii=1_idef,INT(minmax_buffer_size,idef)
minmax = aval_params%min_p - REAL(ii,rp)*(aval_params%max_p- &
aval_params%min_p)/100_rp
if (minmax.GT.0.0_rp) then
min_p = minmax
end if
end do
! buffer at maximum p boundary
max_p = aval_params%max_p + minmax_buffer_size*(aval_params%max_p- &
aval_params%min_p)/100_rp
! buffer at minimum pitch angle boundary
if (aval_params%min_pitch_angle.GE.korc_zero) then
do ii=1_idef,INT(minmax_buffer_size,idef)
minmax = aval_params%min_pitch_angle - REAL(ii,rp)* &
(aval_params%max_pitch_angle-aval_params%min_pitch_angle)/100_rp
if (minmax.GT.0.0_rp) then
min_pitch_angle = minmax
end if
end do
else
min_pitch_angle = aval_params%min_pitch_angle
end if
! buffer at maximum pitch angle boundary
do ii=1_idef,INT(minmax_buffer_size,idef)
minmax = aval_params%max_pitch_angle + REAL(ii,rp)* &
(aval_params%max_pitch_angle-aval_params%min_pitch_angle)/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(P_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
call RANDOM_NUMBER(rand_unif)
P_buffer = min_p + (max_p - min_p)*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,aval_params%dR)
Z_test = Z_buffer + random_norm(0.0_rp,aval_params%dZ)
P_test = P_buffer + random_norm(0.0_rp,aval_params%dp)
T_test = T_buffer + random_norm(0.0_rp,aval_params%dth)
! Test that pitch angle and momentum are within chosen boundary
do while ((T_test .GT. aval_params%max_pitch_angle).OR. &
(T_test .LT. aval_params%min_pitch_angle))
T_test = T_buffer + random_norm(0.0_rp,aval_params%dth)
end do
do while ((P_test.LT.aval_params%min_p).OR. &
(P_test.GT.aval_params%max_p))
P_test = P_buffer + random_norm(0.0_rp,aval_params%dp)
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)
! Repetedly put location of sample into first entry of spp%vars%Y,
! beacause spp%vas%X isn't filled until the end of this subrouting.
! This allows us to use the machinery for calling PSPLINES
! interpolation routines here also.
spp%vars%Y(1,1)=R_buffer
spp%vars%Y(1,2)=0
spp%vars%Y(1,3)=Z_buffer
call get_profiles(params,spp%vars,P,F)
! write(output_unit_write,'("ne",E17.10)') spp%vars%ne(1)
! write(output_unit_write,'("Te",E17.10)') spp%vars%Te(1)
! write(output_unit_write,'("Zeff",E17.10)') spp%vars%Zeff(1)
! Update avalanche parameters with interpolated fields to be used
! in call to avalanche distribution function
call update_avalanche_params(params,spp%vars)
fRE0=fRE(cos(deg2rad(T_buffer)),P_buffer)
else
psi0=psi1
fRE0=fRE1
end if
psi1=PSI_ROT(R_test,spp%Ro,spp%sigmaR,Z_test,spp%Zo, &
spp%sigmaZ,theta_rad)
spp%vars%Y(1,1)=R_test
spp%vars%Y(1,2)=0
spp%vars%Y(1,3)=Z_test
call get_profiles(params,spp%vars,P,F)
call update_avalanche_params(params,spp%vars)
fRE1=fRE(COS(deg2rad(T_test)),P_test)
! write(output_unit_write,'("psi0: ",E17.10)') psi0
! write(output_unit_write,'("psi1: ",E17.10)') psi1
! 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)* &
sin(deg2rad(T_test))*fRE1/(R_buffer*EXP(-psi0)* &
sin(deg2rad(T_buffer))*fRE0)
if (ratio .GE. 1.0_rp) then
R_buffer = R_test
Z_buffer = Z_test
P_buffer = P_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
P_buffer = P_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,aval_params%dR)
Z_test = Z_buffer + random_norm(0.0_rp,aval_params%dZ)
P_test = P_buffer + random_norm(0.0_rp,aval_params%dp)
T_test = T_buffer + random_norm(0.0_rp,aval_params%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,aval_params%dth)
end do
do while ((P_test.LT.min_p).OR.(P_test.GT.max_p))
if (P_test.lt.0) then
P_test=abs(P_test)
exit
end if
P_test = P_buffer + random_norm(0.0_rp,aval_params%dp)
end do
psi0=psi1
fRE0=fRE1
psi1=PSI_ROT(R_test,spp%Ro,spp%sigmaR,Z_test,spp%Zo, &
spp%sigmaZ,theta_rad)
spp%vars%Y(1,1)=R_test
spp%vars%Y(1,2)=0
spp%vars%Y(1,3)=Z_test
call get_profiles(params,spp%vars,P,F)
call update_avalanche_params(params,spp%vars)
fRE1=fRE(COS(deg2rad(T_test)),P_test)
ratio = indicator(psi1,psi_max_buff)*R_test*EXP(-psi1)* &
sin(deg2rad(T_test))*fRE1/(R_buffer*EXP(-psi0)* &
sin(deg2rad(T_buffer))*fRE0)
if (ratio .GE. 1.0_rp) then
R_buffer = R_test
Z_buffer = Z_test
P_buffer = P_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
P_buffer = P_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. &
(P_buffer.LE.aval_params%max_p).AND. &
(P_buffer.GE.aval_params%min_p).AND. &
(T_buffer.LE.aval_params%max_pitch_angle).AND. &
(T_buffer.GE.aval_params%min_pitch_angle)) THEN
R_samples(ii) = R_buffer
Z_samples(ii) = Z_buffer
P_samples(ii) = P_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(P_samples,spp%ppp,MPI_REAL8, &
mom,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
spp%vars%g = SQRT(1.0_rp + mom**2)
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'
DEALLOCATE(mom)
if (params%mpi_params%rank.EQ.0_idef) then
DEALLOCATE(R_samples)
DEALLOCATE(Z_samples)
DEALLOCATE(PHI_samples)
DEALLOCATE(P_samples)
DEALLOCATE(T_samples)
end if
end subroutine Avalanche_4D