Kinetic particle simulations in a global toroidal geometry

TL;DR

The paper presents an upgraded gyrokinetic code for global toroidal simulations, incorporating kinetic particle pushers and validating ion orbit loss and zonal flow phenomena in the DIII-D tokamak.

Contribution

It introduces a fully kinetic particle pusher and coupling of core and edge regions in a global gyrokinetic simulation framework.

Findings

Ion orbit loss fractions increase with ion energy.

Simulation results agree with theoretical and experimental observations.

Successful verification of the new code's capability in tokamak edge and core phenomena.

Abstract

The gyrokinetic toroidal code (GTC) has been upgraded for global simulations by coupling the core and scrape-off layer (SOL) regions across the separatrix with field-aligned particle-grid interpolations. A fully kinetic particle pusher for high frequency waves (ion cyclotron frequency and beyond) and a guiding center pusher for low frequency waves have been implemented using cylindrical coordinates in a global toroidal geometry. The two integrators correctly capture the particle orbits and agree well with each other, conserving energy and canonical angular momentum. As a verification and application of this new capability, ion guiding center simulations have been carried out to study ion orbit losses at the edge of the DIII-D tokamak for single null magnetic separatrix discharges. The ion loss conditions are examined as a function of the pitch angle for cases without and with a radial electric field. The simulations show good agreement with past theoretical results and with experimentally observed feature in which high energy ions flow out along the ion drift orbits and then hit the divertor plates. A measure of the ion direct orbit loss fraction shows that the loss fraction increases with the ion energy for DIII-D in the initial velocity space. Finally, as a further verification of the capability of the new code, self-consistent simulations of zonal flows in the core region of the DIII-D tokamak were carried out. All DIII-D simulations were perfomed in the absence of turbulence.

Tables

Table 1. Table 1: Initial conditions for FK integrator

Parameter	Core	Cross-Separatrix
$R/R_0$	1.154	1.289
$Z/R_0$	0.0	0.0
$\zeta$	1.570	1.570
$v^R/{\omega }_c{R}_0$	3.371$\times {10}^{-3}$	3.371$\times {10}^{-3}$
$v^Z/{\omega }_c{R}_0$	5.371$\times {10}^{-3}$	5.371$\times {10}^{-3}$
$v^{\zeta }/{\omega }_c$	1.271$\times {10}^{-3}$	1.271$\times {10}^{-3}$
$R_0$ is the $R$ at magnetic axis

Table 2. Table 2: Initial conditions for GC integrator

Parameter	Core	Cross-Separatrix
$R/R_0$	1.147	1.281
$Z/R_0$	0.0	0.0
$\zeta$	1.570	1.570
$v_{\parallel }/{\omega }_c{R}_0$	2.547$\times {10}^{-3}$	3.017$\times {10}^{-3}$
$\sqrt{\mu B_0}/{\omega }_c{R}_0$	4.545$\times {10}^{-3}$	4.642$\times {10}^{-3}$
$B_0$ is the $B$ at magnetic axis

Table 3. Table 3: Parameters for ion orbit loss region

Parameter	Case-I	Case-II
$(R_a/R_0,Z_a/R_0)$	(1.2748, 0.0)	(1.2748, 0.0)
$(R_{in}/R_0,Z_{in}/R_0)$	(0.6134, 0.0)	(0.6134, 0.0)
$(R_{XD}/R_0,Z_{XD}/R_0)$	(0.8442, -0.6680)	(0.6253, -0.6680)
$({\psi }_{XD}/B_0{R}_0^2,{\psi }_a/B_0{R}_0^2)$	(5.253$\times {10}^{-2}$, 4.688$\times {10}^{-2}$)	(5.138$\times {10}^{-2}$, 4.688$\times {10}^{-2}$)

Equations

$$\vec{B}=\nabla \psi (R,Z)\times \nabla \zeta +\frac{F(\psi )}{R}\hat{\zeta },$$

$$B_R=-\frac{1}{R}\frac{\partial \psi }{\partial Z},B_Z=\frac{1}{R}\frac{\partial \psi }{\partial R},B_{\zeta }=\frac{F(\psi )}{R},$$

$$J^{-1}=\nabla R\cdot (\nabla \zeta \times \nabla Z).$$

$$x=R\cos \zeta ,$$

$$y=R\sin \zeta ,$$

$$z=Z.$$

$$\vec{v}=v^R{\vec{e}}_R+v^{\zeta }{\vec{e}}_{\zeta }+v^Z{\vec{e}}_Z,$$

$$\vec{\mathbb{E}}=-\nabla\phi=-\biggl{[}\frac{\partial\phi}{\partial R}\nabla R+\frac{\partial\phi}{\partial\zeta}\nabla\zeta+\frac{\partial\phi}{\partial Z}\nabla Z\biggr{]},$$

$$v^R=\dot{R},v^{\zeta }=\dot{\zeta },v^Z=\dot{Z},$$

$${\vec{e}}_R=\cos \zeta \hat{x}+\sin \zeta \hat{y},$$

$${\vec{e}}_{\zeta }=-R\sin \zeta \hat{x}+R\cos \zeta \hat{y},$$

$${\vec{e}}_Z=\hat{z}.$$

$$f(\psi )=f(1,i)+f(2,i)\Delta \psi +f(3,i)\Delta {\psi }^2,$$

$$\Delta\psi=\frac{\psi_{Lim}}{N_{\psi}-1}\Rightarrow N_{sep}=\texttt{integer}\biggl{[}\frac{\psi_{sep}}{\Delta\psi}\biggr{]}\Rightarrow\Delta\psi=\frac{\psi_{sep}}{N_{sep}},$$

$$f(R,Z)=f(1,i,j)+f(2,i,j)\Delta R+f(3,i,j)\Delta R^2+$$

$$f(4,i,j)\Delta Z+f(5,i,j)\Delta R\Delta Z+f(6,i,j)\Delta Z\Delta R^2+$$

$$f(7,i,j)\Delta Z^2+f(8,i,j)\Delta R\Delta Z^2+f(9,i,j)\Delta R^2\Delta Z^2,$$

$$\displaystyle\Delta r=\frac{1}{{N_{SOL}-1}}\biggl{\{}R_{LF}(\psi=\psi_{1},Z=0)-$$

$$\displaystyle R_{LF}(\psi=\psi_{sep},Z=0)\biggr{\}},$$

$$\displaystyle\psi_{SOL}(i_{SOL})=\psi\biggl{(}R_{LF}(\psi=\psi_{sep},Z=0)+$$

$$\displaystyle(i_{SOL}-1)\times\Delta r,Z=0\biggr{)},$$

$$N_{core}=\texttt{integer}\biggl{[}\frac{R_{LF}(\psi=\psi_{sep},Z=0)-R_{0}}{\Delta r}\biggr{]}$$

$$\psi_{0}=\psi\biggl{(}R_{LF}(\psi=\psi_{sep},Z=0)-(N_{core}-1)\times\Delta r,Z=0\biggr{)},$$

$$\displaystyle\psi_{core}(i_{core})=\psi\biggl{(}R_{LF}(\psi=\psi_{0},Z=0)+$$

$$\displaystyle(i_{core}-1)\times\Delta r,Z=0\biggr{)},$$

$$\displaystyle N_{private}=\texttt{Integer}\biggl{[}\frac{R_{LF}(\psi=\psi_{sep},Z=Z_{0})}{\Delta r}-$$

$$\displaystyle\frac{R_{LF}(\psi=\psi_{pmax},Z=Z_{0})}{\Delta r}\biggr{]}$$

$$\displaystyle\psi_{p0}=\psi\biggl{(}R_{LF}(\psi=\psi_{sep},Z=0)-$$

$$\displaystyle(N_{private}-1)\times\Delta r,Z=Z_{0}\biggr{)},$$

$$\displaystyle\psi_{private}(i_{private})=\psi\biggl{(}R_{LF}(\psi=\psi_{p0},Z=Z_{0})+$$

$$\displaystyle(i_{private}-1)\times\Delta r,Z=Z_{0}\biggr{)},$$

$$\frac{dR}{dS}=b_R,$$

$$\frac{dZ}{dS}=b_Z,$$

$$N_i=\text{integer}(\frac{l_i}{\Delta S_{i0}}),$$

$$\Delta S_i=\frac{l_i}{N_i-1}.$$

$$\biggl{[}\frac{\partial}{\partial t}+\vec{v}\cdot\nabla+\frac{q_{c}}{m}(\vec{\mathbb{E}}+\vec{v}\times\vec{B})\cdot\frac{\partial}{\partial\vec{v}}\biggr{]}f_{FK}=0,$$

$$\frac{d}{dt}\vec{r}=\vec{v},\quad\quad\frac{d}{dt}\vec{v}=\frac{q_{c}}{m}\biggl{[}\vec{\mathbb{E}}+\vec{v}\times\vec{B}\biggr{]}.$$

Full text

I Introduction

One of the important challenges in achieving a viable operating regime for ITER and future fusion reactors is associated with the nonlinear turbulent dynamics of the plasma in the scrape-off layer (SOL) loarte2007power . The plasma characteristics in SOL can greatly affect the overall confinement properties of the device and also regulate the heat load to the tokamak wall. It can also influence the level of fusion ash, the impurity dynamics, sheath physics and plasma shaping effects. Furthermore, the SOL dynamics can degrade the current drive performance of radio frequency (RF) waves, through their impact on the density threshold conditions for the onset of parametric decay instabilities cesario2010current . An in-depth undertanding of the mechanisms determining the width of the SOL layer remains an outstanding open problem. A study of the SOL plasma dynamics is challenging due to the multiple spatial and temporal scales associated with different energy sources (instabilities) in that region. Fluid simulation transport codes such as UEDGE wising1996simulation and SOLPS rozhansky2009new are normally used to simulate the SOL dynamics. These fluid codes use a set of fluid transport equations that are based on the Braginskii equations braginskii1965review . However, the results show a number of discrepancies between experimental findings and fluid simulations, especially in the characteristics of the radial electric field, parallel ion flow, impurity radiation, etc. chankin2007discrepancy ; erents2004comparison ; wischmeier2011assessment . It is believed that kinetic effects could be a significant contributor in the SOL to processes like ion orbit losses pan2014co , X point losses chang2002x , nonlocal turbulent transport omotani2013non , plasma sheath dynamicswesson1995effect , parametric decay instabilities cesario2010current ; kuley2009stabilization ; kuley2010parametric ; kuley2010lower , etc. To correctly model many of these effects, one requires a kinetic approach that covers the closed and open field line regions across the separatrix, and includes the realistic SOL physics and tokamak geometry. Due to the difficulty in the accessibility of diagnostics in the SOL region, such global kinetic simulations can help develop useful insights for predicting the plasma dynamics in that region for present and future reactors such as ITER and DEMO. A laudable effort in this direction has been the development of the massively parallel kinetic simulation code XGC-1 chang2009compressed that takes an approach based on first-principles and has emerged as an efficient method for describing the complex physics of turbulent transport. Another widely used and successful tool, the gyrokinetic toroidal code (GTC) has undergone continuous development for the past two decades and has been applied to the study of plasma transport in the core region lin1998turbulent . GTC is a well-benchmarked, first-principles code which has been extensively applied to the investigation of neoclassical transport lin1997neoclassical ; dong2016effects , microturbulence xie2017new ; fulton2016gyrokinetic ; xiao2009turbulent ; holod2016effect ; xiao2015gyrokinetic ; zhang2012nonlinear , mesoscale Alfv$\acute{e}$n eigenmodes zhang2012nonlinear ; spong2012verification ; wang2013radial excited by energetic particles, macroscopic MHD modes mcclenaghan2014verification ; liu2014verification ; bao2017conservative (kink and tearing modes) and radio frequency (RF) waves kuley2013verification ; bao2014particle ; kuley2015verification ; kuley2015nonlinear ; bao2015global ; bao2016electromagnetic ; bao2016nonlinear in the core region. However, the assumptions used in studying turbulence in the core region may not be valid in the SOL region. GTC normally uses conventional magnetic flux coordinates, in which the equations of motion encounter a mathematical singularity of the metric on the magnetic separatrix surface. This is due to vanishing of the poloidal magnetic field near X-point as well as singular behaviors of safety factor and Jacobian of Boozer coordinates near separatrix. Recently, the GTC was extended to separately study instabilities in the SOL and core regions of a field reversed configuration (FRC) using Boozer coordinates. However, the code still did not have the capability to couple these two regionsfulton2016gyrokinetic ; fulton2016bgyrokinetic . The difficulty lay in the discontinuity of the poloidal angle across the separatrix, in the Boozer coordinates. This limitation restricted the code’s usage to electrostatic simulations either in the core or the SOL region, with no cross-separatrix coupling.

In our present work we report a significant enhancement of the GTC code, called global toroidal code using X point (GTC-X) through the development of a new global nonlinear particle simulation model that couples the tokamak core and SOL regions. The code also provides a realistic treatment of the separatrix region through the Equilibrium Fitting (EFIT) lao1990equilibrium ; ren2011high and IPREQSrinivasan equilibrium data files generated by experimental discharges. A particular feature of GTC-X is the use of a cylindrical coordinate system for the advancement of the particle dynamics, which allows particle motion in arbitrary shaped flux surfaces including the magnetic separatrix and the magnetic X-point in the tokamak. Currently, GTC-X has both fully kinetic and guiding center particle dynamics, but XGC-1 has only guiding center particle dynamics. The long-term goal of GTC-X is to develop electromagnetic (EM) global simulation model for coupling core and SOL using both guiding center and fully kinetic particle dynamics. GTC-X in its present form is an electrostatic code, which has been applied to simulate ion temperature gradient instability in the FRC geometry, which has no toroidal magnetic field Bao19 .

As a first step in developing this nonlinear particle simulation model, we have developed a method of field-aligned particle-grid interpolations using an axisymmetric mesh in the cylindrical coordinates, which takes advantage of the smallest number of grid points in the direction of the magnetic field with high resolution in any given poloidal plane. These field-aligned particle-grid interpolations can achieve the same numerical efficiency as the field-aligned mesh in magnetic flux cooordinates, employed earlier in the global code GTC xiao2015gyrokinetic , and in the flux-tube codes such as FENICA hariri2013flux and GEM Sturdevant2017 . The gain in computational efficiency by using appropriate coordinates and computational mesh helps to optimize turbulence simulations of large devices like ITER and DEMO. We have also developed a fully kinetic (FK) particle pusher to capture the effect of high frequency waves (ion cyclotron frequency and beyond) and a guiding center (GC) pusher to describe the particle dynamics associated with low frequency waves (much smaller than the ion cyclotron frequency).

To test the effectiveness of these enhancements and appropriately benchmark the code, we have carried out ion guiding center simulations to study ion orbit losses at the edge of the DIII-D tokamak for single null magnetic separatrix discharges. Using model calculations, some analytic expressions of such losses have been presented by Miyamoto miyamoto1996direct and have been used in the past to estimate the loss region in velocity space for JET, JT-60 and ITER. Stacey has introduced the effect of ion orbit loss and X point loss in his fluid calculations for the interpretation of fluid transport in the edge region stacey2013effect ; stacey2011effect . We apply the new simulation model and benchmark GTC-X against some of these past results. We have examined the ion orbit losses as a function of the pitch angle, both in the presence and absence of an electric field. The simulations show good agreement with past theoretical results and with many experimentally observed features. A measure of the ion direct orbit loss fraction shows that the loss fraction increases with the ion energy for DIII-D in the initial velocity space.

As a testimony to the capability of GTC-X to reproduce the existing physical phenomena, we carried out self-consistent simulations of zonal flows in the core region of the DIII-D tokamak. The collisionless damping of the zonal electric field to a non-zero steady-state value, verifies the famous theory of Rosenbluth and Hinton Rosenbluth on the collisionless damping of zonal flows.

The paper is organized as follows. Section II contains a detailed description of the coordinate system and the representation of the equilibrium and fluctuating field quantities in the code. Section III describes the computational mesh used in the simulations. The physics model explaining the dynamics of the particles is described in section IV followed by our simulation results for the ion orbit losses in section V. The verification of the damping of zonal flows is given in Section VI. Section VII provides a brief summary and some concluding discussions.

II Representation of equilibrium and coordinate system

In describing the collective dynamics of a tokamak plasma, it is convenient to represent the physical variables in terms of equilibrium and fluctuating quantities. The equilibrium is described by the Grad-Shafranov equation. The fluctuations representing wave activity contribute to plasma transport. The input equilibrium magnetic field configuration for GTC-X can be generated in any of the equilibrium solvers, EFIT or IPREQ, and is normally expressed as a function of the magnetic flux function. In the MHD approximation (ignoring guiding center orbit effects) the plasma profiles of temperature, density, current etc. are also functions of the flux function. In the GTC-X version we use cylindrical coordinates $(R,\zeta,Z)$ where $R$ is the distance from the geometry axis, $\zeta$ is the toroidal angle and $Z$ is the direction of the torus symmetry axis. The representation of the magnetic field for an axisymmetric system is then given by,

[TABLE]

where $\psi(R,Z)$ is the poloidal flux function that labels the magnetic surfaces for both closed and open field lines [cf. Fig 1(b)] and $F(\psi)$ is the poloidal current function [cf. Fig 1(a)], which provides the components of magnetic field in cylindrical coordinates as:

[TABLE]

The Jacobian for this system can be written as

[TABLE]

The cylindrical toroidal coordinate system is related to the standard Cartesian system as follows

[TABLE]

By defining contravariant basis vectors $\vec{e}^{R}=\nabla R$, $\vec{e}^{\zeta}=\nabla\zeta$, $\vec{e}^{Z}=\nabla Z$ the velocity and the electric field can be written as

[TABLE]

[TABLE]

where

[TABLE]

The equilibrium inputs from EFIT only provide equilibrium quantities on a coarse mesh, which usually contains a few tens of grid points in the $R$ and $Z$ directions. However, the micro-scale turbulence demands much denser grid points in $R$ and $Z$ directions. Therefore, it is necessary to map the equilibrium mesh to a dense computational mesh, in order to achieve sufficient numerical accuracy. For a 1D function $f(\psi)$, such as the poloidal current function $F(\psi)$, we can use the following B-spline representation xiao2015gyrokinetic

[TABLE]

where $f(2,i)$ and $f(3,i)$ are coefficients related to the first and second order differential in $\psi$ direction, which are calculated using finite difference method on a spline mesh. In our simulations we have calculated the grid size in $\psi$ using the following three steps:

[TABLE]

where $\psi_{Lim}$ is the poloidal flux function at the limiter point, $N_{\psi}$ is the number of grid points in $\psi$, $N_{sep}$ is the grid point number at the separatix, and $\psi_{sep}$ is the poloidal flux function at the separatrix. This calculation will provide an accurate value of $\psi$ at the separatrix in the simulation. However, it will give some minor differences in calculating the $\psi$ at the limiter points.

A 2D function $f(R,Z)=\sum_{n}\mathbb{F}_{n}(R)\mathbb{G}_{n}(Z)$ can be expressed as

[TABLE]

where $\Delta R=R_{i+1}-R_{i}$, $i=1,2,3,\cdots,nR$, and $\Delta Z=Z_{j+1}-Z_{j}$, $j=1,2,3,\cdots,nZ$. Eq. (10) is derived by using the 1D B-spline functions of $\mathbb{F}_{n}(R)=\mathbb{F}_{n}(1,i)+\mathbb{F}_{n}(2,i)\Delta R+\mathbb{F}_{n}(3,i)\Delta R^{2}$ and $\mathbb{G}_{n}(Z)=\mathbb{G}_{n}(1,j)+\mathbb{G}_{n}(2,j)\Delta Z+\mathbb{G}_{n}(3,j)\Delta Z^{2}$. The spline coefficients $f(1:9,i,j)$ are calculated from the spline coefficients $\mathbb{F}_{n}(1:3,i)$ and $\mathbb{G}_{n}(1:3,j)$. Fig. 1(a) and Fig. 1(b) represent the equilibrium poloidal current function on a uniform flux grid and poloidal flux function on the rectangular R-Z grid respectively, for DIII-D shot $\#$158103 at 3050 ms PhysRevLett.114.105002 . In GTC-X we use these two functions in Eq.(2) to calculate the magnetic field components for DIII-D [cf. Fig. 2]. In the particle pusher, we interpolate these field quantities at the particle position using a 2D spline function, as described in Eq. (10). The magnetic field is divergence-less if we use the expressions of $B_{R},B_{Z}$, and $B_{\zeta}$ as given by Eq.(2). The equilibrium poloidal flux $\psi$ is numerically calculated by a 2D spline representations. Hence for this axisymmetric system, the $\nabla\cdot\vec{B}=0$ is guaranteed numerically in the simulation.

III Field aligned mesh for fluctuating quantities

The field aligned mesh has maximum numerical efficiency and accuracy to address the nonlinear physics of drift wave turbulence. Because the particle moves much faster in the direction parallel to the magnetic field than when it drifts across the magnetic field, the parallel wave vector is usually much smaller than the perpendicular wave vector for the drift wave instabilities. Thus, one only requires a small number of grid points to resolve the parallel wave length, which greatly saves the computational costs and suppresses the numerical high $k_{||}$ modes efficiently. Furthermore, it helps to simplify the implementation of the field solver, since the field aligned mesh can exactly decouple the parallel and perpendicular directions.

III.1 Radial grid:

GTC code has in the past utilized the field aligned mesh in Boozer coordinates and has successfully applied it to gyrokinetic simulations of drift waves in the core region of toroidal plasmas lin1998turbulent . In this work, we extend the original GTC field aligned mesh from Boozer coordinates in the core region to cylindrical coordinates in the whole domain across the separatrix. In order to create the field aligned mesh for the whole tokamak domain, we build equi-spaced radial grids (spacing of size $\Delta r$) on the outer mid-plane for each field line, which can be calculated as:

[TABLE]

where $R_{LF}$ is the radial position as a function of poloidal flux and $Z$ on the low field side, and $N_{SOL}$ is the number of field lines in the SOL region. $\psi_{1}$ is the maximum value of the poloidal flux on the limiter points (plasma facing components). GTC-X also has the capability to use a nonuniform grid, with the grid size in the perpendicular direction correlated with the local gyro-radius. The poloidal flux at each field line of the simulation mesh in the SOL region is:

[TABLE]

where $i_{SOL}\in[1,N_{SOL}]$.

Then the field line number in the core region $N_{core}$ can be calculated as:

[TABLE]

where $R_{0}$ is the major radius. The poloidal flux at the innermost grids in the core region $\psi_{0}$ can be determined as:

[TABLE]

where $\psi_{0}$ is not necessarily equal to the poloidal flux value at the magnetic axis. The poloidal flux of each field line for the simulation mesh in the core region is:

[TABLE]

where $i_{core}\in[1,N_{core}]$.

Similarly, the field line number in the private region (i.e., the space/dome below separatrix X-point, which has no connection of flux lines with main plasma) $N_{private}$ is calculated as:

[TABLE]

where $\psi_{pmax}$ is the poloidal flux at the innermost flux surface in the private region, and $Z_{0}$ is the minimum value in $Z$ direction. The poloidal flux at the innermost grids in the private region $\psi_{p0}$ can be determined as:

[TABLE]

where $\psi_{p0}$ is not necessarily equal to $\psi_{pmax}$. The poloidal flux at each field line of the simulation mesh in the private region is:

[TABLE]

where $i_{private}\in[1,N_{private}]$.

III.2 Poloidal grid:

After this, we trace each field line along the poloidal direction and calculate the length of the poloidal projection of each field line by using the unit magnetic vector as follows:

[TABLE]

and

[TABLE]

Here, $b_{R}=-(1/B_{p}R)(\partial\psi/\partial Z)$, $b_{Z}=(1/B_{p}R)(\partial\psi/\partial R)$ and $B_{p}$ is the poloidal magnetic field. The starting points for tracing the field lines are on the low field side of outer mid-plane ($Z=0$ plane) for the core region and, on the low field side of the $Z=Z_{0}$ plane for the SOL and the private regions. The step $\Delta S_{t}$ for tracing each field line is at least $10^{-5}$ smaller than the order of magnitude of the length of the poloidally projected field line, and we correct the new position for each tracing step by using Newton’s method, with the error of the order of floating precision, in order to ensure that it lies on the same flux surface,. All the advanced positions during the tracing of the poloidal field line are recorded as tracing grids. The outermost simulation boundary is determined by the limiter, as shown by the brown line. As a result, the tracing grids outside the limiter in the SOL and private region are removed by using ray casting algorithm wiki.org . In the closed field line region, we calculate the length of the poloidal projection of the field line for each flux surface. However, in the open field line regions, we set the intersection points of the field lines and the limiter as the simulation boundary points and then calculate the length of the poloidal projection of the field line between each of the two intersection points on a single continuous field line.

The number of simulation grids for each continuous poloidal field line is calculated as,

[TABLE]

where $l_{i}$ is the poloidal field line length of i-th field line and $\Delta S_{i0}$ is the approximate grid size along the poloidal direction. The exact poloidal grid size $\Delta S_{i}$ is then determined by $l_{i}$ and $N_{i}$ as:

[TABLE]

The positions of the field aligned grid along poloidal direction can be derived from poloidal grid size $\Delta S_{i}$, the much smaller tracing step $\Delta S_{t}$ and the positions of the dense tracing grids.

In the cylindrical coordinate representation, we only create the field aligned mesh in the poloidal direction, and the meshes on different poloidal plane are identical (as shown in Fig. 3) due to axisymmetry. In order to decouple the exact parallel and perpendicular directions and use a small number of poloidal planes, the particle-grid gathering and scattering operations can be done exactly along the parallel direction by interpolating the value on the identical field line.

IV Physics model for particle dynamics

The efficiency of particle simulation strongly depends on the advancement of the dynamical quantities. In the GTC-X we have developed particle pushers for both fully kinetic particles as well as guiding center particles using cylindrical coordinates $(R,\zeta,Z)$ in a global toroidal geometry. The models for the fully kinetic dynamics and guiding center dynamics and, the numerical methods associated with the time advancement of the physical quantities (particle position and guiding center) are described in the following sections.

IV.1 Fully kinetic particle dynamics

Fully kinetic particle dynamics is described by the six dimensional Vlasov equation

[TABLE]

where $f_{FK}$ is the fully kinetic particle distribution function, $q_{c}$ is the particle charge, and $m$ is the particle mass. The above equation governs the evolution of $f_{FK}$ in the fully kinetic approach.

The time evolution of the phase space coordinates of a single particle, in the presence of a self-consistent electromagnetic field, is governed by the Lorentz-force equation as follows:

[TABLE]

In our simulation we compute the marker particle trajectory Eq. (24) using the time centered Boris push method kuley2013verification ; kuley2015verification ; Boris ; Birdsall ; Tajima ; wei2015method as discussed in section 4.2. The Lagrangian for the single particle motion in the cylindrical coordinates is written as

[TABLE]

where $\vec{A}$ is the magnetic vector potential. Now the components of generalized momenta are

[TABLE]

In our simulation we use the poloidal flux function $\psi$, rather than the vector potential $\vec{A}$ to represent the equilibrium magnetic field components for calculating the toroidal canonical angular momentum. Two constants of motion for fully kinetic single-particle dynamics in cylindrical coordinates are defined by Sturdevant2017 :

• •

Kinetic Energy

[TABLE]

• •

Toroidal angular momentum from Eq. (26) and Eq. (1)

[TABLE]

since $RA_{\zeta}=\psi$

IV.2 Boris push for fully kinetic particle dynamics

Boris scheme is the most widely used orbit integrator in explicit particle-in-cell (PIC) simulation of plasmas. In this paper we have extended our Boris push scheme in GTC from Boozer coordinates kuley2015verification to cylindrical coordinates. This scheme offers second order accuracy while requiring only one force (or field) evaluation per step. The interplay between the PIC cycle and the Boris scheme is schematically represented in Fig.2 of Ref.kuley2015verification . At the beginning of each cycle the position of the particles and their time centered velocity $\vec{v}(t-1/2)$ as well as the grid based electromagnetic fields $\vec{\mathbb{E}}(t),\vec{B}(t)$ are given.

At the first step, we add the first half of the electric field to update the velocity from $(t-1/2)$ to $t$ as follows:

[TABLE]

One may write the components of velocity at particle position at $t$ as

[TABLE]

where $\alpha=R,\zeta,Z$. For an orthogonal cylindrical system the above Eq. (28) can be rewritten as as follows:

[TABLE]

where $A_{1}=\cos(\zeta_{2}-\zeta_{1}),$ $B_{1}=\cos\zeta_{1}\sin\zeta_{2}-\sin\zeta_{1}\cos\zeta_{2}$, $\zeta_{1}=\zeta(t-1/2)$, $\zeta_{2}=\zeta(t)$, $R(1)=R(t-1/2)$ and $R(2)=R(t)$.

In the second step we consider the rotation of the velocity at time $(t)$. Rotated vector can be written as

[TABLE]

where $\vec{T}=(q_{c}\vec{B}/m)(\Delta t/2)$ and $\vec{s}=2\vec{T}/(1+T^{2})$. The components of rotated vector become

[TABLE]

[TABLE]

[TABLE]

where $Q=(q_{c}/m)(\Delta t/2)$, $P=2/(1+T^{2})(q_{c}/m)(\Delta t/2)$, $g_{\zeta\zeta}=R^{2}$, and $J^{2}=det(\texttt{g}_{\alpha\beta})=R^{2}$. In the third step, we add the other half electric acceleration to the rotated vectors to obtain the velocity at time $(t+1/2)$

[TABLE]

To update the particle position we need to recover $\vec{v}(t+1/2)$, which can be done through the following transformation (cf. Fig.2 of Ref.kuley2015verification dark purple arrow)

[TABLE]

where $\gamma=R,\zeta,Z$. However, the basis vector $\nabla\gamma(t+1/2)$ is still unknown, since $\gamma(t+1/2)$ does not exist in standard leap-frog scheme. Here we use a prediction for $\gamma(t+1/2)$ as

[TABLE]

After we find the velocity at time $(t+1/2)$, we can update the particle position using the leap-frog scheme as

[TABLE]

In Eq. (35) we have the dot-product of two basis vectors at different time steps. We have calculated this equation in the similar fashion as discussed in Eqs. (28)- (29).

In this section, we have described the time advancement of the dynamical quantities such as velocity and position of the particle using time centered approach. However, the self-consistent simulation requires update of particle weight (representing perturbed distribution function), guiding center and electric field. In our future simulations, we will use the second order Runge-Kutta (RK) method, to advance these quantities kuley2015verification ; kuley2015nonlinear .

IV.3 Guiding center dynamics

Guiding center particle dynamics is described by the five-dimensional phase space

[TABLE]

where $f_{GC}$ is the guiding center distribution function, $\vec{X}$ is the guiding center position and $v_{\parallel}$ is the parallel velocity. The evolution of the guiding center distribution function can be described by the following equations of guiding center motionbrizard2007foundations :

[TABLE]

where $\vec{B}^{*}=\vec{B}+Bv_{\parallel}/\omega_{c}\nabla\times\hat{b}$, $\mu=mv_{\perp}^{2}/2B$, and $B_{\parallel}^{*}=\hat{b}\cdot\vec{B}^{*}$. The $\vec{\mathbb{E}}\times\vec{B}$ drift velocity $\vec{v}_{\mathbb{E}}$, the grad-B drift velocity $\vec{v}_{g}$ and curvature drift velocity $\vec{v}_{c}$ are given by

[TABLE]

In GC description following order is adopted:

[TABLE]

where $\omega$ is the frequency of the mode of interest, $k_{\parallel}$ and $k_{\perp}$ are the wave vectors in the parallel and perpendicular direction, respectively.Two constants of motion for guiding center dynamics are defined by:

• •

Total Energy

[TABLE]

• •

Toroidal angular momentum

[TABLE]

For an axisymmetric system, Eq. (39) can be rewritten in cylindrical coordinate $(R,\zeta,Z)$ as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

For guiding center particle dynamics GTC normally uses second order Runge Kutta (RK) method (cf. Fig.2 of Ref. kuley2015verification ).

To test our integration schemes (FK and GC), we consider the collisionless single particle motion of the trapped one in the core and cross-separatrix regions using the initial conditions from Table 1 and Table 2 for DIII-D geometry, respectively. The projection of the FK and GC trapped particle orbit on R-Z plane in the core and cross-separatrix regions are shown in Fig. 3. Both integrators correctly capture the trapped particle orbit and agree well with each other. In Fig. 4, the temporal variation of the toroidal angle $\zeta$ of these orbits is shown. The bounce motion of the particles is implied by the oscillations of $\zeta$ and the gradual change in the mean level of these oscillations with time is an evidence of the toroidal precession of these orbits. The conservation properties of our integrators are tested with two exact constant of motion, viz., kinetic energy E and toroidal angular momentum $p_{\zeta}$.

In numerical simulation the dynamical quantities accumulate error with time. Fig. 5 shows the relative variation of E and $p_{\zeta}$ for FK Boris integrator over a time of 15000 cyclotron period for different time step sizes in core and cross-separatrix regions. Boris algorithm maintains adequate accuracy with 20 time steps per cyclotron period $(\omega_{c}\Delta t=0.292)$. The relative energy error arises mostly due to floating point cutoff (floating precision). Whereas, the GC second order RK demonstrate that E and $p_{\zeta}$ can converge with 240 time steps per bounce period $(\omega_{b}\Delta t=0.027)$ in core and cross-separatrix regions, where $\omega_{b}$ is the bounce frequency of the particle [cf. Fig.( 6)]. From these convergence studies, it is found that the FK Boris integrator provides better energy convergence than the second order Runge Kutta GC integrator. If there are high frequency electromagnetic perturbations, the condition of $\omega_{per}\Delta t<1$ will set an upper bound for the time step size. In the present simulations, the particles orbit are studied without electric and perturbed magnetic field. However, previously Wei et al.wei2015method have demonstrated the effect of electric field on trapped particle orbit (Ware pinch) in the core region of tokamak using Boozer coordinates.

As the particles move around, some of them happen to leave the domain of simulation. This, after a certain duration of time, would result in erroneous ensemble averages. To improve the accuracy of our simulations, we brought the escaping particles back into the simulation domain using the following procedure: We first found the flux at the position of the particle at a given instant of time $\psi_{p}(t)$, the simulation domain being the interval $(\psi_{1},\psi_{2})$, and the poloidal angle at that instant $\theta_{p}(t)$. If $\psi_{p}(t)<\psi_{1}$ or $\psi_{p}(t)>\psi_{2}$ i.e. if the particle lay outside the simulation domain, we re-initialized $\psi_{p}$ with its value at the earlier time step when the particle was within the simulation domain (i.e. $\psi_{p}(t)=\psi_{p}(t-\Delta t)$) and $\theta_{p}$ was re-initialized as $\theta_{p}(t)=2\pi-\theta_{p}(t-\Delta t)$. This technique preserves the nature of the motion of trapped particles at the domain boundaries. Knowing $\psi_{p}(t)$ and $\theta_{p}(t)$, we obtained the new coordinates $(R_{p},Z_{p})$ of the guiding center by finding expressions for $R$ and $Z$ in terms of $\psi$ and $\theta$. This was done by constructing 2D splines for $R=R(\psi,\theta)$ and $Z=Z(\psi,\theta)$ in the same way as mentioned in Section II.

[TABLE]

[TABLE]

Here $\Delta\psi=\psi_{i+1}-\psi_{i}$ and $\Delta\theta=\theta_{j+1}-\theta_{j}$ where $i=1,2,...,N_{1}$ and $j=1,2,...,N_{2}$. $N_{1}$ and $N_{2}$ depend on the density of the spline grids that our accuracy of interpolation demands. The velocity of the escaping particle is re-initialized with its value at the previous time-step when the particle lay within the simulation domain. This was done to ensure the conservation of the total energy and the toroidal angular momentum. In order to test the above technique, we applied it to the cross-separatrix particle in Fig.  3, using the separatrix as the boundary of our simulation domain such that the particle was brought back every time it crossed the separatrix. The relative error in the total energy and the canonical angular momentum for different time-steps is shown in Fig.  7. We observe that the conservation of the invariants is better at the smaller time-step sizes ($\omega_{b}\Delta t=0.014,0.027$) than that at a bigger time-step size ($\omega_{b}\Delta t=0.055$). However, the overall extent of conservation is not as good as observed earlier, in the absence of a domain boundary. This is because the boundary condition used here is not a perfect one. It however gives us good results in the self-consistent simulations of large number of particles over a long time intervals, as we will observe later in Section VI, where we verify zonal flows in the core region.

V Ion orbit loss near plasma edge

In a tokamak plasma, ions have drift motions due to the gradient-B and curvature drift. As a result, ion orbits are shifted from a magnetic surface. Due to this shift of ion orbit from magnetic surface, the hot ions that exist close to the separatrix, can pass near the X point region. In this region, the poloidal magnetic field is very weak and the ions have a very small poloidal displacement in time. These ions experience vertical curvature and grad-B drifts and move towards the divertor, resulting in ion orbit loss pan2014co ; miyamoto1996direct . In an axisymmetric configuration, the drift orbit of ion is obtained by solving the following three equations as described in Sec 4.3. These are the conservation of kinetic energy, magnetic moment and canonical angular momentum.

[TABLE]

[TABLE]

[TABLE]

The subscript ‘a’ means the value at the starting point of the ions. The conservation of above quantities with $B_{\zeta}\simeq B$ gives the drift orbit surface as follows miyamoto1996direct :

[TABLE]

[TABLE]

where $\xi_{0}=(v_{\parallel a}/v_{a})$ is the initial pitch angle, and $\xi_{0}>0$ denotes the ion moving in the co-current direction. In the above equation, we assume that $B_{\zeta}>0$, and use the approximation $B_{\zeta}\propto(1/R)$. A minus sign must be used instead of a plus sign in Eq. (49) when $B_{\zeta}<0$. We also consider that the poloidal magnetic field is directed clockwise and the ion toroidal drift $(\vec{v}_{g})$ is approximately vertically up and away from the X point. The coordinates of the X point are $(R_{X},Z_{X})$ and the poloidal magnetic flux function of X point satisfies $(\partial\psi/\partial R)_{X}=0$ and $(\partial\psi/\partial Z)_{X}=0$. When the initial position $R_{a}$ of a test ion is at the outer midplane, that is $R_{a}>R_{X}$ there are following two cases for orbit ion loss [cf. Fig. 8(a)-(b)] miyamoto1996direct :

When $\xi_{0}>0$ and $R_{a}>R_{X}$, the relative position of magnetic surface and drift orbit surface [cf. Eq. (50)], clearly indicates that direct ion orbit loss is not possible. However, when the initial position of the test ion is selected as $R_{a}<R_{X}$ the following Fig. 8(c) is possible for ion orbit loss.

In the above GTC-X simulations (cf. Fig. 8) we did not include the radial electric field on ion orbit loss. However, one must mention here that in a tokamak scenario the radial electric field shifts the velocity space boundaries separating trapped and passing orbits. Also the ion orbit loss contributes to the generation of radial electric field, which affects the ion orbit loss. Therefore, the ion orbit loss and the radial electric field should be considered self-consistently for more accurate calculation. Similarly, it is important to study the effect of ion orbit loss and X point loss for understanding the plasma transport stacey2011effect , particle distribution stacey2013effect , plasma rotation pan2014co , etc. Finally, the ion orbit loss region in our calculation is based on the assumptions that the ions are collisionless and the loss region is not filled by the collisions or turbulence. Simulation related to these effects will be reported in a future work. To see how the electric field affects the velocity space boundaries, we have calculated the minimum velocity of the ion and particle loss fractions for DIII-D in the following paragraph.

In order to find the minimum initial ion speed $v_{0}$ that is required for the ion to reach the final location (near X point) we combine Eq. (50) and Eq. (51). The minimum loss speed of ion is

[TABLE]

where $\Psi=(q_{c}/m)(\psi_{a}-\psi_{XD})$ and $\tilde{\phi}=q_{c}(\phi_{a}-\phi_{XD})/(m\omega_{c}^{2}R_{0}^{2}/2)$. For passing ion orbit loss the sign of the second term of the numerator in Eq. (52) is positive and for the trapped ion orbit loss it is negative, respectively. We solve Eq. (52) to find the ion orbit loss region for DIII-D model parameters as described in Table 3. Here, we examine the loss of ion for which $-1\leq\xi_{0}\leq 0$. The loss region in the initial velocity space of an ion starting from outer midplane point $R_{a}$ is described by the solid curves in Fig. 9(a). Case-I and Case-II represent boundary curves CD (C*′D′) and AB (A′B′) of the loss region without $\tilde{\phi}$ (with $\tilde{\phi}$), respectively. Boundary curves BC (B′C′*) i.e, the condition for mirror reflection of the test ion without $\tilde{\phi}$ (with $\tilde{\phi}$) is given by Case-I(ii),where $\tilde{\phi}>0$ and $\tilde{\phi}<0$ correspond to the outward and inward radial electric field, respectively. The value of minimum energy $(mv_{0}^{2}/2)$ is obtained from Eq. (52) [cf. Fig. 9(a)] without $\tilde{\phi}$ demonstrate good agreements with the simulation results of Fig. 8(a) and Fig. 8(b).

After the minimum loss energy is determined, particle loss fraction due to direct ion orbit loss is calculated. The ion loss rate is determined by the rate of supply of ions to the loss regions. The corresponding cumulative particle loss fraction of ions, which follow distribution function $g$ in velocity space is defined as miyamoto1996direct :

[TABLE]

For a Maxwellian distribution function the above equation turns out to be, $f_{L}=\int_{-1}^{1}\Gamma(3/2,\epsilon_{min}(\xi_{0}))d\xi_{0}/2\Gamma(3/2)$, where $\epsilon_{min}(\xi_{0})=mv_{0}^{2}(\xi_{0})/2T_{i}$ is the energy corresponding to the minimum velocity for which ion orbit loss is possible [cf. Eq. (52)], $T_{i}$ is the ion temperature, and $\Gamma(3/2,\epsilon_{min}(\xi_{0}))$ is the upper incomplete gamma function of order 3/2. The $\xi_{0}$ dependent cumulative particle loss fraction $f_{L}$ starting from $R_{a}>R_{X}$ without $\tilde{\phi}$ (ABCD) and with $\tilde{\phi}$ (A′B′C′D′) are shown in Fig. 9(b). These results clearly demonstrate that radial electric field shifts the velocity space boundaries of loss regions for trapped and passing particles. Fig. 9(c) represents the cumulative ion loss fraction for different ion energy ($T_{i}$). Pan et al., have carried out detailed studies of minimum loss energy and cumulative loss fraction for the ions starting from different poloidal and radial positions pan2014co . However, for verification and application of our new capability we have carried out the above calculation for an ion starting from the outer mid-plane only.

VI Verification of Zonal flows in the core region

In this section, we demonstrate the collisionless damping of zonal flows and verify the Rosenbluth-Hinton theory Rosenbluth by observing the variation of the zonal electric field $E(\psi)$, at a given point in the core region, with respect to time. Zonal flows are low frequency electrostatic modes that are spontaneously generated by turbulence and, in turn, play an important role in regulating the turbulence. They are important in the study of transport processes in the core region. In order to simulate these zonal flows, we solve the flux-surface averaged gyrokinetic Poisson’s equation xiao2015gyrokinetic at every time step.

[TABLE]

Here $\langle...\rangle$ implies flux-surface average. $T_{i},n_{i}$ and $Z_{i}$ are the temperature, density and atomic number of the ion species respectively. $e$ is the unit electronic charge and $n_{e}$ is the electron density. $\rho_{i}$ is the ion gyroradius. Since we are dealing with zonal quantities, we are only concerned with the component of the Laplcian perpendicular to the flux surfaces. As shown in Xiao et al.xiao2015gyrokinetic , the flux-surface averaged Laplacian can be written in the flux coordinates $(\psi,\theta,\zeta)$ as

[TABLE]

where $J_{0}(\psi)$ is the flux-surface averaged Jacobian. In writing Eq.(55) we have neglected the higher order toroidal coupling terms by assuming a small inverse-aspect-ratio. In our simulations we load a uniform ion temperature and density profile. As a result, $T_{i}$ and $n_{i}$ are constant all throughout the simulation domain. Integrating Eq.(54) once with respect to $\psi$, we obtain the zonal electric field as follows:

[TABLE]

Here $\langle\rho_{c}\rangle=e(Z_{i}\langle\overline{n}_{i}\rangle-\langle n_{e}\rangle)$ is the flux-averaged total charge density. $\overline{n}_{i}$ is the ion guiding center density. Solving Eq.(56) at every time step, we have incorporated the field information in our GC pusher which updates the particle position accordingly. However, Eq.(56) gives the zonal electric as a function of $\psi$ while the guiding center equations of motion in cylindrical coordinates (Eqs. 41-44) require the $R$ and $Z$ components of the zonal field. Hence we simply do a coordinate transformation to obtain $E_{R}$ and $E_{Z}$ i.e. $E_{R}=(\partial\psi/\partial R)E_{\psi}$ and $E_{Z}=(\partial\psi/\partial Z)E_{\psi}$, where $E_{\psi}=E(\psi)$.

In our simulation of zonal flows, we started off by loading the particles (here, thermal ions) in a certain domain within the core region such that, they are uniformly distributed between any two consecutive flux surfaces. For this, we had to calculate the volume associated with each flux surface. We began by taking an approximation of the volume associated with each grid cell and then added them up over the respective flux surfaces to get the flux surface volumes. After that we loaded the particles accordingly and calculated the number of particles lying within a certain neighborhood of each flux surface (i.e. their weights associated with a given flux surface), thereby obtaining the flux surface volumes in terms of the number of particles. The newly calculated flux surface volumes would also serve as the flux-averaged Jacobian $J_{0}(\psi)$ at each flux surface. As the system evolves with time, the particles leaving the simulation domain are brought back using the reflective boundary conditions discussed earlier, in order to maintain quasi-neutrality and energy conservation. Full-f ions and adiabatic electrons are used in the present study of zonal flows. The initial velocities of the particles followed a Maxwellian distribution. The simulation domain extended from $\psi_{1}=0.4\psi_{sep}$ to $\psi_{2}=0.9\psi_{sep}$ and the domain-width along the outer midplane was $\Delta\approx 68.02\rho_{i}$.

The zonal flows were initiated by driving the system for a certain duration of time with a low magnitude electric field having a sinusoidal dependence on $\psi$. After switching off the external perturbation (i.e. the driving field) the system was allowed to evolve self-consistently with time. Fig. 10 shows the time evolution of the zonal field $E_{\psi}$ after removing the external perturbation. Finally $E_{\psi}$ settles down to a residual value, in the absence of collisions. This non-zero value is the Rosenbluth-Hinton residual level and is of great importance in determining the level of turbulence in the plasma Rosenbluth .

VII Summary and Conclusion

As a first step in developing the self-consistent global simulation model to couple the SOL and the core regions in tokamaks by incorporating the sepratrix, we have developed the particle dynamics for FK and GC particles using cylindrical coordinates. To get the maximum numerical efficiency, field aligned mesh in the cylindrical coordinates is developed for the whole device from magnetic axis to the material wall (plasma facing components). These field-aligned particle-grid interpolations using an axisymmetric mesh in cylindrical coordinates helps to avoid the difficulty associated with the X point. Finally, we have extended this particle orbit simulation techniques to study the ion orbit loss near the X-point of DIII-D single null divertor. Based on three constants of motion, the minimum loss speed and ion orbit loss fraction at the edge of the tokamak are calculated. However, these particles need proper boundary conditions near the divertor wall. Finally, we also verified the physics of zonal flows using the equlibrium of the DIII-D tokomak. Presently GTC-X does not have the capability of handling the plasma particle loss to the material wall (plasma facing components), recycling the neutral particles, and calculating the Monte-Carlo neutral particle transport with a charge exchange and ionization interaction with the plasma.

VIII Acknowledgements

A.K. would like to thank Dr. R. Nazikian and Ms. D. Sharma for providing EFIT and IPREQ data. This work is supported by Board of Research in Nuclear Sciences (BRNS Sanctioned $\#$ 39/14/05/2018-BRNS), Science and Engineering Research Board EMEQ program (SERB Sanctioned $\#$ EEQ/2017/000164), US Department of Energy (DOE) SciDAC ISEP Program, China National Magnetic Confinement Fusion Energy Research Program, (Grant No. 2013GB111000) and IISc startup grant. A.S. thanks the Indian National Science Academy (INSA) for their support under the INSA Senior Scientist Fellowship scheme. All the simulations were carried out on Cray XC40 - SahasraT at Supercomputing Education and Research Centre (SERC), IISc.