@brief Subroutine that generates an exponentially decaying radial distribution in an elliptic torus as the initial spatial condition of a given particle species in the simulation. @details As a first step, we generate an exponentially decaying radial distribution in a circular cross-section torus as in \ref korc_spatial_distribution.exponential_torus. Then we transform this spatial distribution to a one in an torus with an elliptic cross section, this following the same approach as in \ref korc_spatial_distribution.elliptic_torus.
@param[in] params Core KORC simulation parameters. @param[in,out] spp An instance of the derived type SPECIES containing all the parameters and simulation variables of the different species in the simulation. @param fl Variable used in the Newton-Raphson method for finding the radial position of each particle. @param fr Variable used in the Newton-Raphson method for finding the radial position of each particle. @param fm Variable used in the Newton-Raphson method for finding the radial position of each particle. @param rl Variable used in the Newton-Raphson method for finding the radial position of each particle. @param rr Variable used in the Newton-Raphson method for finding the radial position of each particle. @param rm Variable used in the Newton-Raphson method for finding the radial position of each particle. @param relerr Tolerance used to determine when to stop iterating the Newton-Raphson method for finding . @param rotation_angle This is the angle in \ref korc_spatial_distribution.elliptic_torus. @param r Radial position of the particles . @param theta Uniform deviates in the range representing the uniform poloidal angle distribution of the particles. @param zeta Uniform deviates in the range representing the uniform toroidal angle distribution of the particles. @param X Auxiliary vector used in the coordinate transformations. @param Y Auxiliary vector used in the coordinate transformations. @param X1 Auxiliary vector used in the coordinate transformations. @param Y1 Auxiliary vector used in the coordinate transformations. @param pp Particle iterator.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(KORC_PARAMS), | intent(in) | :: | params | |||
| type(SPECIES), | intent(inout) | :: | spp |
subroutine exponential_elliptic_torus(params,spp)
TYPE(KORC_PARAMS), INTENT(IN) :: params
TYPE(SPECIES), INTENT(INOUT) :: spp
REAL(rp) :: fl
REAL(rp) :: fr
REAL(rp) :: fm
REAL(rp) :: rl
REAL(rp) :: rr
REAL(rp) :: rm
REAL(rp) :: relerr
REAL(rp), DIMENSION(:), ALLOCATABLE :: rotation_angle
REAL(rp), DIMENSION(:), ALLOCATABLE :: r
REAL(rp), DIMENSION(:), ALLOCATABLE :: theta
REAL(rp), DIMENSION(:), ALLOCATABLE :: zeta
REAL(rp), DIMENSION(:), ALLOCATABLE :: X
REAL(rp), DIMENSION(:), ALLOCATABLE :: Y
REAL(rp), DIMENSION(:), ALLOCATABLE :: X1
REAL(rp), DIMENSION(:), ALLOCATABLE :: Y1
INTEGER :: pp
ALLOCATE(X1(spp%ppp))
ALLOCATE(Y1(spp%ppp))
ALLOCATE(X(spp%ppp))
ALLOCATE(Y(spp%ppp))
ALLOCATE( rotation_angle(spp%ppp) )
ALLOCATE( theta(spp%ppp) )
ALLOCATE( zeta(spp%ppp) )
ALLOCATE( r(spp%ppp) )
! Initial condition of uniformly distributed particles on a
! disk in the xz-plane
! A unique velocity direction
call init_u_random(10986546_8)
call init_random_seed()
call RANDOM_NUMBER(theta)
theta = 2.0_rp*C_PI*theta
call init_random_seed()
call RANDOM_NUMBER(zeta)
zeta = 2.0_rp*C_PI*zeta
! Uniform distribution on a disk at a fixed azimuthal theta
call init_random_seed()
call RANDOM_NUMBER(r)
do pp=1_idef,spp%ppp
rl = 0.0_rp
rr = spp%r_outter
fl = fzero(rl,spp%r_outter,spp%falloff_rate,r(pp))
fr = fzero(rr,spp%r_outter,spp%falloff_rate,r(pp))
if (fl.GT.korc_zero) then
relerr = 100*ABS(fl - fr)/fl
else
relerr = 100*ABS(fl - fr)/fr
end if
do while(relerr.GT.1.0_rp)
rm = 0.5_rp*(rr - rl) + rl
fm = fzero(rm,spp%r_outter,spp%falloff_rate,r(pp))
if (SIGN(1.0_rp,fm).EQ.SIGN(1.0_rp,fr)) then
rr = rm
else
rl = rm
end if
fl = fzero(rl,spp%r_outter,spp%falloff_rate,r(pp))
fr = fzero(rr,spp%r_outter,spp%falloff_rate,r(pp))
if (fl.GT.korc_zero) then
relerr = 100*ABS(fl - fr)/fl
else
relerr = 100*ABS(fl - fr)/fr
end if
end do
r(pp) = rm
end do
Y = r*SIN(theta)
X = r*COS(theta) + spp%shear_factor*Y
rotation_angle = 0.5_rp*C_PI - ATAN(1.0_rp,1.0_rp + spp%shear_factor);
X1 = X*COS(rotation_angle) - Y*SIN(rotation_angle) + spp%Ro
Y1 = X*SIN(rotation_angle) + Y*COS(rotation_angle) + spp%Zo
spp%vars%X(:,1) = X1*SIN(zeta)
spp%vars%X(:,2) = X1*COS(zeta)
spp%vars%X(:,3) = Y1
DEALLOCATE(X1)
DEALLOCATE(Y1)
DEALLOCATE(X)
DEALLOCATE(Y)
DEALLOCATE(rotation_angle)
DEALLOCATE(theta)
DEALLOCATE(zeta)
DEALLOCATE(r)
end subroutine exponential_elliptic_torus