Subroutine for generating a uniform elliptic torus as the initial spatial condition of a given particle species in the simulation.
An initial spatial distribution following the uniform distribution of torus is modified through a shear transformation and a rotation to generate a uniform spatial distribution on tori with elliptic cross sections. First, we obtain the uniform spatial distribution in a torus of minor radius , see torus. Then, we perform a shear transformation that changes the cross section of the torus from circular to a tilted ellipse. In cylindrical coordinates this shear transformation is given by:
where is the shear factor of the transformation. Here, and are the radial and vertical position of the particles uniformly distributed in a circular torus, and are their new positions when following a uniform distribution in a torus with elliptic circular cross section. The center of the ellipse is , and , where and is the center of the initial circular torus. The major and minor semi-axes of the tilted ellipse cross section is:
Finally, we rotate the ellipse cross section anticlockwise along by , so the major semi-axis is parallel to the -axis.
Modify this approximation.
| 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 elliptic_torus(params,spp)
!! @note Subroutine for generating a uniform elliptic torus as the initial
!! spatial condition of a given particle species in the simulation. @endnote
!! An initial spatial distribution following the uniform distribution of
!! [[torus]] is modified through a shear transformation and a rotation to
!! generate a uniform spatial distribution on tori with elliptic cross sections.
!! First, we obtain the uniform spatial distribution in a torus of minor radius
!! \(r_0\), see [[torus]]. Then, we perform a shear transformation that changes
!! the cross section of the torus from circular to a tilted ellipse. In
!! cylindrical coordinates this shear transformation is given by:
!!
!! $$R' = R + \alpha Z,$$
!! $$Z' = Z,$$
!!
!! where \(\alpha\) is the shear factor of the transformation.
!! Here, \(R\) and \(Z\) are the radial and vertical position of the particles
!! uniformly distributed in a circular torus, \(R'\) and \(Z'\) are their
!! new positions when following a uniform distribution in a torus with
!! elliptic circular cross section. The center of the ellipse is
!! \(R_0' = R_0 + \alpha Z_0\), and \(Z_0 = Z_0\), where \(R_0\) and \(Z_0\)
!! is the center of the initial circular torus. The major and minor semi-axes
!! of the tilted ellipse cross section is:
!!
!! $$a' = \left[ - \frac{2r_0^2}{\alpha \sqrt{\alpha^2 + 4} - (2+\alpha^2)}
!! \right]^{1/2},$$
!! $$b' = \left[ \frac{2r_0^2}{\alpha \sqrt{\alpha^2 + 4} + (2+\alpha^2)}
!! \right]^{1/2}.$$
!!
!! Finally, we rotate the ellipse cross section anticlockwise along
!! \((R_0',Z_0')\) by \(\Theta = \cot^{-1}(\alpha/2)/2\), so the major semi-axis is
!! parallel to the \(Z\)-axis.
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 :: rotation_angle
!! This is the angle \(\Theta\) in the equations above.
REAL(rp), DIMENSION(:), ALLOCATABLE :: theta
!! Uniform deviates in the range \([0,2\pi]\) representing the uniform
!! poloidal angle \(\theta\) distribution of the particles.
REAL(rp), DIMENSION(:), ALLOCATABLE :: r
!! Radial position of the particles \(r\).
REAL(rp), DIMENSION(:), ALLOCATABLE :: zeta
!! Uniform deviates in the range \([0,2\pi]\) representing
!! the uniform toroidal angle \(\zeta\) distribution of the particles.
REAL(rp), DIMENSION(:), ALLOCATABLE :: X
!! Auxiliary vector used in the coordinate transformations.
REAL(rp), DIMENSION(:), ALLOCATABLE :: Y
!! Auxiliary vector used in the coordinate transformations.
REAL(rp), DIMENSION(:), ALLOCATABLE :: X1
!! Auxiliary vector used in the coordinate transformations.
REAL(rp), DIMENSION(:), ALLOCATABLE :: Y1
!! Auxiliary vector used in the coordinate transformations.
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)
r = SQRT((spp%r_outter**2 - spp%r_inner**2)*r + spp%r_inner**2)
Y = r*SIN(theta)
X = r*COS(theta) + spp%shear_factor*Y
!> @todo Modify this approximation.
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 elliptic_torus