Development of semi-Lagrangian gyrokinetic code for full-f turbulence simulation in general tokamak geometry

Journal of Computational Physics 283 (2015) 518–540

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Development of semi-Lagrangian gyrokinetic code for full-f turbulence simulation in general tokamak geometry Jae-Min Kwon a,∗ , Dokkyun Yi b , Xiangfan Piao c , Philsu Kim d a

WCI Center for Fusion Theory, National Fusion Research Institute, Daejeon 305-333, Republic of Korea Institute of Mathematical Sciences, Ewha W. University, Seoul, 120-750, Republic of Korea Department of Mathematics, Hannam University, Daejeon 306-791, Republic of Korea d Department of Mathematics, Kyungpook National University, Daegu 702-701, Republic of Korea b c

a r t i c l e

i n f o

Article history: Received 17 March 2014 Received in revised form 8 December 2014 Accepted 10 December 2014 Available online 16 December 2014 Keywords: Semi-Lagrangian Gyrokinetics Turbulence Tokamak Simulation

a b s t r a c t In this work, we report the development of a new gyrokinetic code for full-f simulation of electrostatic turbulence in general tokamak geometry. Backward semi-Lagrangian scheme is employed for noise-free simulation of the gyrokinetic Vlasov equation with finite Larmor radius effects. Grid systems and numerical interpolations are implemented to deal with arbitrary equilibrium information given in a GEQDSK format file. In particular, we introduce an adaptive interpolation technique for fluctuating quantities, which have elongated structures along equilibrium magnetic fields. This field-aligned interpolation can reduce the required number of grid points to represent the fluctuating quantities. Also, it is shown that the new interpolation allows us to choose bigger time sizes with better simulation accuracy. Several benchmark simulation results are presented for comparison with previously known cases. It is demonstrated that the new code can reproduce the well known results of zonal flow and linear ITG instabilities in concentric circular equilibrium. It is also shown that the code can capture the effects of plasma shaping on the zonal flow and ITG instabilities, and the stabilization effects of the shaping result in significant up-shifts of the threshold of ion temperature gradient for ITG in nonlinear simulation. © 2014 Elsevier Inc. All rights reserved.

1. Introduction In magnetic fusion machine such as tokamak, various plasma turbulences with relatively low frequencies (i.e. ω  Ωi , where Ωi denotes the ion gyro-frequency) are responsible for the anomalous loss of confined thermal energy and particles. For the improvement of the machine performance, it is crucial to understand the physics of plasma turbulence and numerical simulation is an essential tool for the study. Gyrokinetic modeling of strongly magnetized plasma has been very useful in describing the physics of such low frequency turbulences [1]. Since the early pioneering work on the subject [2], there have been enormous efforts to improve numerical methods for gyrokinetic simulation with better accuracy and more realistic conditions: we refer Refs. [3–8,10], which is certainly not a complete list of all related works. Recently, gyrokinetic code developments have been mainly focused on the capabilities to simulate self-consistently the dynamics of plasma turbulence and background profile evolution in more realistic geometry, in particular without any assumption on scale separation.

*

Corresponding author. E-mail address: [email protected] (J.-M. Kwon).

http://dx.doi.org/10.1016/j.jcp.2014.12.017 0021-9991/© 2014 Elsevier Inc. All rights reserved.

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

519

Assessing these goals, both particle-in-cell (PIC) and continuum approaches have their own merits and issues. The PIC approach provides a simulation method with good conservation of key physical quantities and ideal parallelization scalability, which is an important feature for efficient simulation on massively parallel supercomputers in these days. However, the PIC method has issues of small scale noises, which often contaminate physical results and sometimes even render the results useless. The continuum approach is free from these noise issues and guarantees high quality simulation results, if sufficient spatio-temporal resolutions are provided. Employing well developed numerical algorithms in computational fluid dynamics area, continuum gyrokinetic codes are in both rapid development and wide utilization in the analysis of experimental data in these days. However, numerical dissipations pertaining to the continuum schemes require careful control of the discretization algorithm and resolution of simulation [11], which often leads to the issue of computational costs. There have been efforts to combine the advantages of both approaches, and the so-called particle-continuum [12] or semi-Lagrangian methods [13] were proposed for gyrokinetic simulation. In these approaches, the evolution of gyrokinetic distribution function is traced by integrating the equations of motion for gyrokinetic particles (or equivalently the characteristic equations of gyrokinetic equation) as in the PIC simulations. To resolve the small scale noise issues, the distribution function of the particles is reset periodically on a fixed Eulerian grid system, which then filters out the noises systematically. The semi-Lagrangian schemes can be categorized into two classes according to the time integration direction of the characteristic equations i.e. forward and backward schemes [14]. In the forward semi-Lagrangian (FSL) schemes, the characteristic equations are integrated in forward time direction as the same way in conventional PIC methods, which allows the employment of high order time integration schemes and provides a good flexibility in the time step control of simulation. However, the reset (or reconstruction) of the distribution function poses issues in the discretization of the phase space for gyrokinetic equation. Since the particles are “randomly scattered” in the phase space at a target time for the reset, the spatial interpolation of the distribution function becomes a challenging job, especially if we consider the high dimensionality (D = 5) of the phase space. The backward semi-Lagrangian (BSL) schemes resolve this interpolation issue by finding the values of the distribution function at uniformly spaced grid points. This is achieved by taking the uniform grid points as the initial conditions for the integration of the characteristic equations. By tracing the characteristic equations backward in time, we can find the departure points of the characteristic curves at present time t = tn , which arrive at the uniform grid points at future time step t = tn+1 . Then, using the invariance of the distribution function along the characteristics i.e. df /dt = 0, the values of the distribution function at the future time step can be found for the uniform grid points. Since all interpolations are performed in uniform grid system, many known interpolation techniques can be readily employed. However, it should be noted that this advantage in the spatial interpolation comes with a cost in time integration. Since the initial conditions for the time integration of the characteristics are placed at the future time step, while the distribution function is known at the present time step, the time discretization of the equations leads to a highly nonlinear set of equations involving self-consistent electrostatic field satisfying the gyrokinetic Poisson equation, and numerical procedures such as fixed point iterations are required to solve the nonlinear equations. In this paper, we report the development of a full-f gyrokinetic code based on the BSL scheme. The new code solves an electrostatic version of gyrokinetic equation [15]. In the code development, we put special emphasis on the capability to deal with general plasma equilibrium with arbitrary shape of poloidal cross section. Previous studies showed that microinstabilities strongly depend on the poloidal shape of plasma equilibrium [20,21]. Dynamics of zonal flow also depends on the equilibrium geometry [22]. Indeed, there are experimental evidences suggesting a strong dependence of turbulent transport on magnetic configuration [23]. Also, recent theoretical works on the generation of spontaneous plasma rotation pointed out that the up–down asymmetry of plasma equilibrium provides a mechanism to break the symmetry of fluctuation spectrum, which is necessary for the generation of the turbulence driven residual stress and resulting intrinsic torque [24,25]. Considering the importance of all these subjects for better understanding and modeling of turbulent plasma transport in tokamak, it would be useful to develop a new gyrokinetic code which can handle arbitrary plasma equilibrium from experimental data. In this work, we choose the standard GEQDSK format [26] for the equilibrium input and implement all necessary numerical components to deal with the format. The new code will enable us to carry out noise-free simulation of turbulent transport and self-consistent profile evolution in general tokamak geometry, and thus help us analyze and understand experimental data more directly. The remainder of this paper is organized as follows. In Section 2, we explain the simulation model and detailed numerical schemes employed in the code development. The parallelization technique used in the code is explained in Section 3 with discussions on the parallelization scalability of the code. In Section 4, we present several numerical simulation results for benchmark test. Finally, summary and conclusion are given Section 5. 2. Simulation model and numerical scheme 2.1. Simulation model The following form of electrostatic gyrokinetic equation is solved.

dv  ∂ f ∂f dX ∂ f + · + = S src + S sink ∂t dt ∂ X dt ∂ v 

(1)

520

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

Here, f denotes the distribution function of ions in the guiding center space, which depends on time t, the guiding center coordinate X, parallel velocity v  , and magnetic moment μ. S src and S sink represent source and sink term, respectively. The equations for the characteristics are given as follows.

dX dt mi

= v

dv  dt

B∗ B ∗

=−

+

bˆ B ∗

 ×

c qi



μ∇ B + c ∇Φα

B∗  B ∗

· μ∇ B + qi ∇Φα

(2)



(3)

B∗ and B ∗ are defined as B∗ = B + qi v  ∇ × bˆ and B ∗ = bˆ · B∗ , where bˆ is the unit vector along the equilibrium magnetic i field B. c , q i , mi denote the speed of light, ion charge and mass, respectively. The bracket ·α represents the average over the gyro-phase α of each particle i.e. the finite Larmor radius (FLR) effects. For a given axisymmetric equilibrium information, all necessary components of the vectors B and B∗ can be evaluated straightforwardly, and the detailed formulas are given in Appendix A. Electron responses are assumed to be adiabatic. For a given mean electron density ne0 and temperature T e , the total electron density is set as mc



ne = ne0 1 +

e (Φ − Φ f )



(4)

Te

according to an electrostatic potential Φ . · f represents the flux surface average defined as an integral over the poloidal θ and toroidal φ angle on a given flux surface

 2π 0

dθ

Φ f ≡  2π 0

 2π

dθ

0

dφ J Φ

 2π 0

(5)

.

dφ J

Here, J denotes the Jacobian of the coordinate system. The self-consistent electrostatic field Φ satisfies the following form of gyrokinetic Poisson equation.

 

ρ2 1 − ∇ 2 + ∇⊥ · 2ti ∇⊥ Φ + 2 Φ − Φ f = σ λDi λDe  σ = 4π qi dXdv f δ(X + ρ − x)

(6) (7)

√

ρti is defined as ρti = v ti Ωi−1 = T i /mi Ωi−1 with ˆ ˆ ion temperature T i , and ∇⊥ = ∇ − b(b · ∇). For typical tokamak plasma parameters, we have ρti  λDi and the left hand Here, λi and λe denote the ion and electron Debye length, respectively. side of Eq. (8) can be approximated as

 

ρ2 1 − ∇⊥ · 2ti ∇⊥ Φ + 2 Φ − Φ f = σ . λDi λDe

(8)

Eq. (8) is a 3-dimensional (3D) integro-differential equation. Due to the presence of the integral operator · f , application of standard finite difference methods to the equation yield 3D non-sparse matrices, which are hard to handle and solve numerically. However, Eq. (8) can be rewritten as two coupled equations for Φ f and Φ˜ ≡ Φ − Φ f for easier numerical treatments as follows. In an axisymmetric equilibrium, the Jacobian J in the integral is independent of the toroidal angle φ . Then, the flux surface average can be rewritten as

 2π

dθ J Φ¯

Φ f = 0 2π 0

(9)

,

dθ J

where Φ¯ is the toroidal average of Φ (i.e. Φ¯ ≡

 2π 0

dφΦ/2π ). Now, by taking the toroidal average of Eq. (8), we have

 

ρ2 1 ¯ f = σ¯ . − ∇⊥ · 2ti ∇⊥ Φ¯ + 2 Φ¯ − Φ λDi λDe

(10)

This is a 2-dimensional (2D) integro-differential equation and much easier to handle than the original 3D equation. If we can solve this 2D equation and obtain Φ¯ , it is straightforward to calculate Φ f . Then, with a known Φ f , we can recast Eq. (8) as a 2D differential equation for Φ˜ = Φ − Φ f as

    ρ2 ρ2 1 − ∇⊥ · 2ti ∇⊥ Φ˜ + 2 Φ˜ = σ + ∇⊥ · 2ti ∇⊥ Φ f . λDi λDe λDi

(11)

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

521

Fig. 1. The radial profiles of source and sink terms.

Since the flux surface average operator is now removed, numerical discretization of this equation will yield highly sparse matrices, which can be efficiently solved by various iterative schemes (detailed numerical schemes to solve Eq. (10) and Eq. (11) are discussed in Appendix B). As a boundary condition for the flux surface averaged electrostatic potential Φ f , we impose Φ f = 0 and Φ f = 0 at

√

the inner χ = χmin and outer χ = χmax radial boundary, respectively. Here, the radial variable χ is defined as χ ≡ ψ/ψ X using the poloidal magnetic flux ψ and the flux value at the last closed flux surface ψ = ψ X . For the remaining component of the potential Φ˜ , we impose a Dirichlet boundary condition as Φ˜ = 0 at the two radial boundaries. The gyro-phase average of the potential can be written as Φα = J 0 (k⊥ ρ )Φk in Fourier space. We follow the procedure developed in Ref. [27] to obtain Φα as follows. Using the Pade approximation of the Bessel function

J 0 (k⊥ ρ ) ≈

1 1 + (k⊥ ρ )2 /4

,

we have the following elliptic problem for Φα with Φ as the source term.



J 0−1 (k⊥ ρ )Φα ≈ 1 −

ρ2 4

 2 Φα = Φ. ∇⊥

(12)

We solve this numerically for each value of μ = ρ 2 Ωi2 mi /2B on the grid points of the perpendicular velocity space for f (the details of the grid system will be explained in the forthcoming paragraphs). The gyro-phase average operator in the source term of the Poisson equation is treated similarly. The sink term in the right hand side (rhs) of Eq. (1) takes the form:

S snk ( f ) = −

A snk (χ )

τsnk



 ∂ 1 ∂2  ( f − f 0) + χ D χ (χ ) ( f − f 0 ) + 2 2 D θ (χ )( f − f 0 ) . χ ∂χ ∂χ χ ∂θ 1 ∂

It is applied to make distribution function f close to a given smooth form f 0 near the radial boundaries to avoid numerical problems. To implement heating source term for flux-driven simulation, we follow the method proposed in Ref. [8]. The source term is numerically expressed as a difference of two Maxwellian distributions f m1 and f m2 with an identical density but different temperature profiles.

S src =

A src (χ )

τsrc

( f m1 − f m2 )

(13)

The size of τsrc and the temperature differences between f m1 and f m2 are adjusted according to the level of input power. Typical shapes of the radial envelopes A snk , D χ , D θ , A src are shown in Fig. 1. 2.2. Numerical scheme The grid system for the gyrokinetic simulation is set based on plasma equilibrium information in the GEQDSK format [26]. A GEQDSK file provides the information of axisymmetric plasma equilibrium, which consists of the distributions of poloidal magnetic flux ψ( R , Z ) and plasma current I ( R , Z ) on a poloidal plane in rectangular coordinate ( R , Z ). Once an equilibrium is given, we can reconstruct the flux surfaces of the equilibrium by finding ψ -contours starting from the last

522

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

Fig. 2. Grid systems for 2D interpolations on poloidal plane.

closed flux surface with ψ = ψ X . During the reconstruction, the inverse mappings R (ψ, θ) and Z (ψ, θ) are also found for ψ and θ = tan−1 ( Z / R ), which will be used to interpolate necessary equilibrium information for the simulation. The spatial grid for the gyrokinetic simulation is generated in the coordinate system (χ , φ, θ), where φ denotes the toroidal angle. The triple (χ , φ, θ) forms a right-handed coordinate system i.e. ∇ χ × ∇φ · ∇θ > 0. The velocity space for gyrokinetic distribution √ √ function is represented by ( v  , u ), where u = 2μ = v ⊥ mi / B. Five-dimensional (5D) uniform grid is generated in the coordinate (χ , φ, θ, v  , u ) for the distribution function. For the calculation of various components of the gyrokinetic equations of motion, which are summarized in Appendix A, we need 2D interpolations of ψ and I , and their first and second derivatives with respect to R , Z . At the start of simulation, reading the equilibrium information, we prepare the interpolations of ψ( R , Z ) and I ( R , Z ) on a uniform ( R , Z )-grid. After finding the flux surfaces and the inverse mappings, another set of interpolations are prepared for R (ψ, θ) and Z (ψ, θ). The local cubic interpolation technique is employed for these interpolations and the details are explained in Appendix B. In the actual implementation, the accuracy and grid resolutions of the interpolations are carefully controlled to ensure that |ψ( R (ψ , θ), Z (ψ , θ)) − ψ |/ψ  1 for all ψ and θ in the simulation domain. Combining these, we can evaluate all the necessary 2D interpolations for the gyrokinetic equations at arbitrary point in the spatial coordinate (χ , θ). Fig. 2 illustrates the schematics of the interpolations in 2D grid systems. For low frequency plasma turbulences described by the gyrokinetic equations, the parallel and perpendicular components of the fluctuations show strong anisotropy k  k⊥ i.e. the fluctuating eddies are highly elongated along B-field. Therefore, the spatial interpolations of the fluctuating potential and gyrokinetic distribution function can be done very efficiently with fewer grid points, if we follow the B-field in finding neighboring poloidal planes. The schematics of the field aligned interpolation is illustrated in Fig. 3, and the detailed procedure is explained in Appendix B. Compared to the conventional interpolation constructed by the tensor product of 1D interpolations in each spatial direction, the field-aligned interpolation can substantially reduce the poloidal and toroidal grid number required to correctly represent the fluctuating potential structures elongated to equilibrium magnetic fields. As will be explained in subsequent paragraphs, the time discretization and integration of the BSL scheme require some numerical procedures to solve a set of highly nonlinear equations. Iteration methods are usually employed to solve these equations, and the reduction of the iteration number is essential for an efficient simulation. The field-aligned interpolation allows us to choose bigger time step sizes and provides better numerical accuracy for the iterations. All these advantages of the field-aligned interpolation are presented and discussed in Appendix B. √ Since the perpendicular velocity variable u = 2μ is an invariant of motion, only the interpolation in v  direction is needed in the velocity space. For the v  interpolation, we employ the local cubic interpolation technique as for the configuration space. We remind that this interpolation plays very crucial roles in determining the detailed physical properties of gyrokinetic distribution function [16,17] e.g. the positivity, conservations of key physical quantities etc. For the temporal discretization of Eq. (1), we employ the Strang operator splitting algorithm. We can write the gyrokinetic Vlasov equation with source and sink terms as

∂f dz = − · ∇ f + S src + S snk ≡ A [ f ] + B [ f ]. ∂t dt

(14)

Here z = (χ , θ, φ, v  , u ), A [ f ] ≡ − dz · ∇ f , and B [ f ] ≡ S src + S snk . If we symbolically denote the time integration of the dt homogeneous part and inhomogeneous part as T A (h) and T B (h) for a time step h, respectively, the time integration over a single time interval t can be written as

T A + B (t ) = T A (t /2) T B (t ) T A (t /2)

(15)

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

523

Fig. 3. Schematics of B-field aligned interpolation on flux surface.

according to the splitting algorithm. If we employ second order accurate time discretizations for both T A and T B , the second order accuracy of T A + B is guaranteed by the splitting. For the inhomogeneous part B [ f ], we apply the Crank–Nicholson algorithm, which is second order accurate. In this work, the homogeneous part A [ f ] is integrated using the conventional second order BSL scheme [18,19]. If we write the characteristic equations of Eq. (14) as dz/dt = U(t , z), the BSL scheme requires to find the departure points of the characteristic curves at present time step t = tn , which arrive at grid points at future time step t = tn+1 = tn + t. Then the values of distribution function at the future time step f n+1 can be interpolated using the present values f n . A second order accurate discretization of the characteristic equations give

z = zn + tU(tn+1/2 , z∗ ), z∗ = zn +

t 2

(16)

U(tn , zn ).

(17)

Here, we have to find zn for a given grid point z and calculate self-consistent electrostatic potential at t = tn+1/2 . These equations are highly nonlinear, and we employ the following form of fixed point iterations to solve them. Step 0. Set initial guess ζ (0) = z − zn , η(0) = z − z∗ , l = 0. Step 1. Update η by calculating η(l+1) = ζ (l) − 2t U(tn , z − ζ (l) ).

Step 2. Set zn+1/2 = z − η(l+1) and interpolate f n+1/2 to calculate Φ(tn+1/2 , z). Step 3. Update ζ by calculating ζ (l+1) = tU(tn+1/2 , z − η(l+1) ). (l+1)

Step 4. Set zn = z − ζ (l+1) and interpolate f n+1 . (l+1) (l) (l) Step 5. Check the convergence of zn by calculating the error ε ≡ |( f n+1 − f n+1 )/ f n+1 |. If the error is within a given tolerance, the iteration is terminated. Otherwise, go back to the Step 1 and increase l by one, and repeat the steps. In the following simulations, the tolerance level is set to ε = 10−6 , which requires 2 ∼ 3 fixed point iterations for the convergence in each time step with t = 3Ω −1 . The initial guesses are taken from the final values of previous time step. 3. Parallelization We employ the MPI library for the parallelization of the new gyrokinetic code. The spatial domain is decomposed in all three spatial directions. Each sub-domain is assigned to one MPI processing element (PE) and interpolations for the distribution function and the characteristic equations are performed on each PE separately. Various information of nearby grid points for the interpolations are shared through MPI communications at all intermediate stages of each time step. To solve the Poisson equation (Eq. (11)) and the elliptic equation for the gyro-phase average (Eq. (12)), we employ the GMRES solver in HYPRE library [29], which is a collection of various parallelized iterative matrix solvers. We implement two interfaces for the solver with different domain decomposition schemes. The default interface for the solver has the same domain decomposition scheme as the other parts of the code i.e. the decomposition in all three spatial directions. The other interface has only radial domain decomposition to reduce the number of PE’s for the GMRES solver and optimize the parallelization efficiency of the solver. In this case, we form a sub-communicator for PE’s belonging to a same flux surface and the source terms for the solver are gathered to the local root PE in the sub-communicator. Then, only those local root PE’s participate in the HYPRE calculation. After solving the equations, each local root PE scatters the solutions to every other PE’s in the same sub-communicator. To analyze the parallelization performance of the new code, we performed two set of scaling tests with a fixed size problem N χ × N θ × N φ × N v  × N u = 192 × 256 × 16 × 32 × 8. Fig. 4 shows the variation of simulation times for increasing

524

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

Fig. 4. Parallelization performances of various parts of the code. Computing times (sec) spent for the interpolations (upper-left), field solvers (upper-right), communications (lower-left), and a single time step (lower-right) are shown. The red curves represent the case with full-domain decomposition for the field solvers. The circular-points curves show the case with radial domain decomposition only for the field solvers.

number of PE’s (N PE ). Here, the interpolation part includes all simulation times spent for the interpolations of 2D equilibrium information, 3D fluctuating field, and 5D distribution function. The field part counts the times consumed by the Poisson and the gyro-phase average solvers. The times for the MPI-communications to exchange nearby grid information are also shown in the figure. The total simulation time includes all of these and additional times spent on diagnosis etc. The first set of simulations (the red curves) were performed with the same domain decomposition algorithm for both the field solvers and other simulation components. As we can see, most computing times are consumed by the interpolations and the solvers. The former part can be efficiently parallelized as we increase N PE as shown in the figure. On the other hand, the time required by the field solvers increases as N PE increases and becomes the most time consuming part as ) are the same with N PE > 1000. In these simulations, the number of PE’s participating in the field solving processes (N PE the total number of PE’s i.e. N PE = N PE . Next, we repeated the tests with the alternative domain decomposition for the field solvers, and the results are shown in the same figure (the green curves). We can see that the time consumed by the communications and the interpolations is limited by the number of radial grid show similar behaviors with the first set of tests. In this set of simulations, N PE = N = 24, 48, 96, 192 for N = 192, 384, 768, 1536, respectively. We can see that the points, and they correspond to N PE χ PE ∼ 200 and increases beyond the number. Clearly, the field solver part is CPU time consumed by the solvers saturates for N PE the main bottleneck of the simulation and further works are needed to improve it and extend the parallelization scalability beyond N PE > 1000. 4. Simulation results In this section, we present several simulation results using two different plasma equilibrium configurations. The first configuration has a concentric circular poloidal cross section with B 0 = 1.9 T, R 0 = 1.3 m, a0 = 0.48 m, as the central magnetic field strength, the major and minor radius, respectively. The q-profile is set as q( ) = 0.84 + 15.99 ×  2 , where  = r / R 0 . The second configuration has a D-shaped poloidal cross section with κ = 1.5, δ = 0.5, s = 0.3 as the elongation, triangularity, and Shafranov-shift, respectively. The second one has the same B 0 = 1.9 T, R 0 = 1.3 m. The radial simulation domain is set as [χmin = 0.2, χmax = 1.0] for both configurations. The poloidal domain is discretized as N χ × N θ = 192 × 256 for the circular case, and N χ × N θ = 192 × 512 (i.e. higher poloidal resolution) for the D-shaped case. In Fig. 5, the poloidal cross sections of the two equilibriums are shown with their q-profiles. First, we performed simulations in axisymmetric configurations with Φ f only, and checked the dynamics of geodesic acoustic mode (GAM) and the level of residual zonal flow (ZF). With flat density and temperature profiles, we put small ini-

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

525

Fig. 5. Poloidal cross section and q-profile of equilibriums used in the simulations.

tial density perturbations to excite Φ f . Fig. 6 shows the temporal behavior of Φ f at χ ≈ 0.6 for the two equilibriums. In the concentric circular case, we can see that the residual ZF level is in good agreement with the Rosenbluth–Hinton formula [30]. Also, the GAM damping rate agrees well with the theoretical prediction as shown in the figure. In the D-shaped case, the GAM damping occurs much faster than the circular case and the residual ZF level is higher than the Rosenbluth–Hinton formula. The refined formula taking into account the plasma shaping effects by Xiao–Catto [22] shows a better match with the simulation result. Next, we performed a set of linear ion temperature gradient (ITG) mode simulations with the local Cyclone parameters: R / L ne = 2.22, R / L T i = 6.91, T e = T i , q = 1.4, sˆ = (r /q)dq/dr = 0.84 at  = 0.18 for the concentric circular equilibrium [31]. In full-f simulation, unlike δ f case, we cannot artificially suppress the radial electric field E r = ∂r Φ f and E r evolves self-consistently depending on the initial choice of distribution function f . If we choose a local Maxwellian as the initial f , a strong E r develops and affects the linear ITG instabilities [33]. Since the Cyclone benchmark results in Ref. [31] were obtained by δ f simulations without any strong E r , we have chosen a canonical Maxwellian [34] as the initial f to avoid the formation of strong E r and its effects on the instabilities. A contour plot of the fluctuating potential Φ˜ is shown in Fig. 7. The growth rate and frequency of the instabilities with different toroidal mode number nφ are also shown in the same figure (note that kθ ρi = nφ qρi /r). In these linear simulations with varying nφ , we choose only one section of equally divided nφ -sections of the full-torus. This one toroidal section is discretized with a given number of toroidal grids N φ = 8. In every time steps, we filter out electrostatic field components with toroidal mode numbers n = 1 in this reduced domain. Then, the remaining n = 1 component provides us the required nφ -component of the potential in the full toroidal domain. Considering the differences in numerical schemes and detailed simulation conditions, the results are in good agreement with Ref. [31]. We also performed a set of linear ITG simulations for the D-shaped equilibrium with varying nφ . As for the circular cases, the initial f is set in a canonical Maxwellian form. Equilibrium profiles for the initial f are chosen to provide key local parameters at  = 0.16 as R / L ne = 2.0, R / L T i = 6.92, T e = T i , q = 1.18, sˆ = 0.91. Fig. 8 shows linear growth rates for different nφ ’s and a typical contour plot of the electrostatic potential with nφ = 14. As we can see, the overall growth rates and frequencies become smaller and the most unstable mode number is shifted toward higher nφ , as compared to the circular case. In axisymmetric toroidal plasmas, the radial electric field E r is determined by the radial force balance equation, in which the poloidal plasma flow enters as a contributor to the determination of E r . The poloidal flow is damped by two different physical mechanisms depending on time scale [32]. For relatively short time scales t ∼ ωb−1 (here ωb represents

526

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

Fig. 6. GAM dynamics and residual zonal flows for the circular (upper) and the D-shaped case (lower).

the bounce frequency of trapped particle), it is mainly determined by the collisionless Landau damping. For long collisional time scales t ∼ νii−1  ωb−1 (here νii denotes the ion–ion collision frequency), it is ultimately determined by the neoclassical poloidal viscosity originating from the Coulomb collisions. The new gyrokinetic code does not include the Coulomb collision operators yet. However, within the short time scales t ∼ ωb−1 , the poloidal flow is mainly damped by the collisionless Landau damping. Then, the radial force balance is quickly established by the neoclassical poloidal viscosity from the Landau damping and a quasi-steady state E r is formed, which can be considered as an approximate equilibrium E r in a collisionless plasma. To investigate the effects of the equilibrium E r on ITG instabilities, we repeated the simulations for the D-shaped equilibrium with a local Maxwellian as the initial f . The same set of the equilibrium profiles for the canonical distribution are used. In this case, strong E r develops and shows oscillatory behaviors in the initial phase of the evolution. To minimize the effects of the transient part of E r on the instabilities, the simulations were proceeded in two steps, as follows. In the first step, with the same initial profiles for ITG simulations, an axisymmetric simulation was performed to find E r and f in a quasi-stationary state (i.e. E r and f were taken at a certain time step after passing through the initially transient states). Then, in the next step, they were used as the initial conditions for ITG simulations with toroidal variations. Fig. 9 shows the growth rates and frequencies of the linear instabilities with different nφ ’s. As a comparison, δ f simulation results using gKPSP code [35] are also plotted in the same figure. In performing the δ f simulations, the same D-shaped equilibrium and initial profiles were used. A fixed E r profile, which is close to the quasi-stationary E r obtained in the first step of the full-f simulation, was imposed during the δ f simulations. Even with the considerable differences in numerical schemes and detailed simulation setup, the full-f and δ f simulation results agree reasonably well with the same trend. We can see that, due to the strong E r , the overall modes are significantly stabilized with higher frequencies, as compared to the cases with the canonical initial distribution function (i.e. without the equilibrium E r ) in the same D-shaped equilibrium (see Figs. 8 and 9). Finally, we performed nonlinear simulations for the two equilibrium cases. The grid resolutions are chosen as N χ × N θ × N φ × N v  × N u = 192 × 256 × 32 × 128 × 8 and 192 × 512 × 32 × 128 × 8 for the circular and the D-shaped case,

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

527

Fig. 7. Linear ITG instabilities in the concentric circular equilibrium. Contour plot of fluctuating potential with kθ ρi = 0.32 (upper) and growth rates and frequencies with different toroidal mode numbers (lower).

respectively. To reduce computational costs, we only simulated a quarter of the toroidal domain (i.e. 0 ≤ φ ≤ π /2). This reduced domain is equally divided by N φ = 32, and then the allowed range of toroidal Fourier modes becomes [−16, 16]. Since we are dealing with only a quarter of the full torus, the effective Fourier mode numbers for the simulations are set as nφ = 0, ±4, ±8, ±12, · · · , ±64. For both equilibriums, the initial distribution function f is set in a local Maxwellian form. The initial ion temperature (T i ) profiles are shown in Fig. 10. The central value of T i is chosen to give a/ρi ≈ 120 at the mid-radius in both cases. We followed the same procedure for the linear ITG simulations in the D-shaped case with local Maxwellian to find quasi-stationary axisymmetric equilibrium f and E r , which are used as the initial conditions for the nonlinear simulations. To prevent too rapid relaxation of the ion temperature profiles, we put 2MW heating source in the core regions as described in Section 2. Fig. 10 shows the initial and final ion temperature and zonal flow profiles for the two equilibrium cases. A temporal sequence of evolving fluctuating potentials Φ˜ is shown for the D-shaped case in Fig. 11. We can see that the quasilinear growth of ITG instabilities is followed by the development of strongly sheared E × B flows in addition to the initial ones, which quench the ITG turbulence largely and lead to the nonlinear saturation. More details of the ion temperature and heat flux evolutions are presented in Fig. 12. Here, the values are averaged over the radial interval χ = [0.45, 0.55]. As we can see, the initially strong heat flux diminishes due to the development of strong E × B flow shear and the ion temperature gradient settles to R / L T i ≈ 5.6 and 6.1 for the circular and the D-shaped case, respectively. We note that these gradient values are significantly higher than the usual linear instability thresholds without E × B flow shears, which can be easily checked by running δ f code. Indeed, linear δ f simulations with fixed nφ = 16 and varying R / L T i in Fig. 13 indicate that the up-shift of the instability threshold by the E × B flow is about 1.5 and 3.1 for the circular and the D-shaped case, respectively. It is noteworthy that the stabilization effect appears more strongly in the D-shaped case, as compared to the circular case. We want to remind that the new gyrokinetic code does not include the Coulomb collisions, which are essential to provide physical dissipations and achieve physically correct steady states of turbulence and transport. Therefore, the states appearing after the turbulence saturation by the axisymmetric shear flows are only quasi-steady, and the profiles evolve slowly. To realize physically correct flux-driven steady states of turbulence and transport, we need to implement the Coulomb collision operators numerically, which will be a main subject of our future studies.

528

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

Fig. 8. Linear ITG instabilities in the D-shaped equilibrium. Contour plot of fluctuating potential with kθ ρi = 0.32 (upper) and growth rates and frequencies with different toroidal mode numbers (lower).

Fig. 9. Linear ITG instabilities for local Maxwellian distribution in the D-shaped equilibrium. For comparison, δ f simulation results with the same equilibrium and profiles are presented.

5. Summary and conclusion In this work, we developed a new gyrokinetic code for full-f simulation of turbulent transport and self-consistent profile evolution in general tokamak geometry. The BSL scheme was employed for the time integration of the electrostatic gyrokinetic equations with the adiabatic electron response model. Implementing various numerical components for the new code,

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

529

Fig. 10. Nonlinear simulation results of the concentric circular (upper) and D-shaped (lower) equilibriums. The initial and final profiles of ion temperature (left) and zonal flow (right) are presented.

Fig. 11. Nonlinear simulation results of the D-shaped equilibrium case. Contour plots of the fluctuating electrostatic potential Φ˜ at different time steps are presented.

we devoted special efforts to deal with general axisymmetric equilibrium given in the GEQDSK format. All 2D interpolations on poloidal plane were performed in a full accordance with the general equilibrium information. The 3D spatial interpolations of distribution function and fluctuating potential were implemented to take into account the property k  k⊥ . We introduced a field-aligned interpolation following ambient equilibrium magnetic field lines. Numerical tests were performed to investigate the advantages of the new interpolation method. It was found that the field-aligned interpolation provides us to choose fewer spatial grid points and bigger time step sizes with better simulation accuracy as compared to the conventional interpolation method. The spatial domain was decomposed in all three directions for the parallelization of the new code. An iterative solver in the HYPRE library was employed to solve the discretized forms of the gyrokinetic Poisson equation and the elliptic equation for the gyro-phase average. To optimize the performance of the solver, an alternative

530

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

Fig. 12. Ion temperature gradients and heat conductivities averaged over by the local gyro-Bohm values.

χ = [0.45, 0.55] from the nonlinear simulations. The conductivities are normalized

Fig. 13. δ f simulation results of the linear ITG threshold of ion temperature gradient. nφ = 16 is chosen.

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

531

domain decomposition was optionally applied. It was found that the performance heavily depends on the number of PE’s working on the solver. Using the concentric circular and D-shaped equilibriums, several benchmark simulations were performed for the comparison with previously known results. The dynamics of axisymmetric potential Φ f were in good agreement with known theories, and the D-shaped case exhibited a higher residual zonal flow level, as compared to the theory for concentric circular geometry. Linear ITG instabilities were calculated with the well known Cyclone parameters, and the results were shown to be in good agreement with the previous benchmark simulation results in the circular case. The results of the D-shaped case indicated that the shaping has a strong stabilization effect on the instabilities. Nonlinear simulations with self-consistent E r showed that the threshold of ion temperature gradient for ITG is significantly up-shifted by the E × B shear suppression effect. The stabilization effect appeared more prominently in the D-shaped case, which indicates the importance of plasma shaping for the confinement improvement in tokamak plasmas. In future works, we will revisit the velocity and configuration space interpolation by employing new numerical techniques and investigate how they affect the various properties of gyrokinetic simulation e.g. the conservation of mass, momentum, and energy. Control of the long term error accumulations in these conserved quantities is also an important topic. As found in this work, there is a close relation between the error from the fixed-point iterations and the interpolation method, which should impact the long term conservation errors. This issue will be further investigated in future works. The Coulomb collision operators will be implemented to study flux-driven turbulence in long term steady states. Also, we will try to improve the parallelization efficiency, especially of the field solver part. The present result indicates that further investigation is needed to find optimal number of PE’s according to the size of simulation, which will be determined by the balance between the reduction of local computations and the increase of the communication costs for the field solvers. Ultimately, we plan to perform transport time scale simulations of plasma turbulences with realistic source models and general equilibrium configuration for the study of various tokamak transport problems. Acknowledgements The authors are grateful to Drs. V. Grandgirard, Y. Sarazin, G. Dif-Pradalier, S. Ku, P. Diamond, X. Garbet, G. Latu for useful discussions. This work was supported by the World Class Institute (WCI) Program of the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT and Future Planning (NRF grant number: WCI 2009-001). The second author was also supported by Priority Research Center Program (NRF grant number: 2009-0093827). The third and last authors were also supported by Basic Science Research Program (NRF grant number 2011-0029013). Appendix A. Guiding center equations of motion in general tokamak geometry In axisymmetric tokamak, an equilibrium magnetic field can be expressed as the following form using the poloidal flux

ψ and plasma current I in the cylindrical coordinate ( R , φ, Z ) 1 ∂ψ

B = ∇φ × ∇ψ + I ∇φ =

R ∂Z

Rˆ +

I R

φˆ −

1 ∂ψ R ∂R

Zˆ ,

(A.1)

and its amplitude is given as



B=

1



R

∂ψ ∂R

2

 +

∂ψ ∂Z

2 + I 2.

(A.2)

Various vector terms involving B in the gyrokinetic equations of motion can be also expressed in terms of ψ and I .

∇ ×B=−

1 ∂I R ∂Z

Rˆ +

1 R



∇2ψ −

1 ∂ψ



R ∂R

φˆ +

1 ∂I R ∂R

Zˆ

(A.3)

∇ × B ˆ ∇B ∇ × bˆ = +b×

(A.4)

bˆ · ∇ × bˆ =

(A.5)

B B·∇ ×B B2

B

    ∂ψ ∂ I 1 1 ∂ψ ∂ψ ∂ I − = 2 2 I ∇2ψ − + R ∂R ∂R ∂R ∂Z ∂Z R B

These expressions can be used to evaluate the terms for the contravariant components of Eq. (2).

∇ψ · B = 0

(A.6)

∇θ · B = ∇θ · ∇φ × ∇ψ ≡ J

(A.7)

∇φ · B =

I

(A.8)

R2

∇ψ · (∇ × B) = −

1 ∂ψ ∂ I R ∂R ∂Z

+

1 ∂ψ ∂ I R ∂Z ∂R

≡

1 R

[ I , ψ] R Z

(A.9)

532

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

∇θ · (∇ × B) = − ∇φ · (∇ × B) = ∇ψ · (∇ × bˆ ) = ∇θ · (∇ × bˆ ) =

1 ∂θ ∂ I

R ∂R ∂Z  1 2

∇ ψ−

R2 1 RB 1

+

1 ∂θ ∂ I R ∂Z ∂R

1 ∂ψ

≡



1 R

[ I , θ] R Z

(A.11)

R ∂R

[ I , ψ] R Z +

I R B2 I

[ψ, B ] R Z

[ I , θ] R Z + [θ, B ] R Z RB R B2   1 1 ∂ψ 1 2

∇φ · (∇ × bˆ ) =

R2 B

(A.10)

∇ ψ−

R ∂R

−

R2 B2

(A.12) (A.13)



∂ψ ∂ B ∂ψ ∂ B + ∂R ∂R ∂Z ∂Z





∂Φ |∇ψ|2 ∂Φ ∇ψ · B × ∇Φ = J I + ∂θ ∂φ R2

∂Φ ∇ψ · ∇θ ∂Φ ∇θ · B × ∇Φ = − J I + ∂ψ ∂φ R2 2

|∇ψ| ∂Φ ∇ψ · ∇θ ∂Φ ∇φ · B × ∇Φ = − − ∂ψ ∂θ R2 R2 ∂ f ∂g

(A.14) (A.15) (A.16) (A.17)

∂g ∂ f

Here, we defined J ≡ ∇ψ × ∇φ · ∇θ and [ f , g ] R Z ≡ ∂ R ∂ Z − ∂ R ∂ Z . Once ψ( R , Z ) and I ( R , Z ) are given in the GEQDSK format, these can be evaluated straightforwardly using numerical interpolations. Appendix B. Grid system and interpolation We discretize the rectangular domain containing the poloidal cross section of plasma equilibrium as

R min = R 0 < R 1 < · · · < R N R = R max ,

R i = R 0 + i R ,

(B.1)

Z min = Z 0 < Z 1 < · · · < Z N Z = Z max ,

Z j = Z 0 + j Z ,

(B.2)

where  R = ( R max − R min )/ N R and  Z = ( Z max − Z min )/ N Z are the uniform spatial mesh sizes in R and Z direction, respectively. The values of ψi , j and I i , j at each grid point are read from a given GEQDSK file. A local cubic interpolation is used to calculate necessary 2D equilibrium information. For a given point ( R , Z ), we find the cell [ R i −1 , R i ] × [ Z j −1 , Z j ] containing the point. Then, we can calculate, for example ψ( R , Z ) as

ψ( R , Z ) =

j +1 

i +1 

c i −i (λ)c j − j (ω)ψi , j ,

(B.3)

i = i −2 j = j −2

where λ = ( R − R i )/ R and

c −2 (x) = c −1 (x) = c 0 (x) = c 1 (x) =

ω = ( Z − Z j )/ Z , and the interpolation coefficients c i, j ’s are given as follows.

x − x3

(B.4)

6 x3 + x2 − 2x 2

2 + x − 2x2 − x3 2 x3 + 3x2 + 2x 6

(B.5) (B.6) (B.7)

We can check the validity and accuracy of the 2D grid systems and interpolations for long term gyrokinetic simulation by performing test particle orbit integration. The gyrokinetic equations of motion in Eqs. (2) and (3) can be integrated for test particles in a shaped axisymmetric equilibrium geometry without the fluctuating potential parts. Then, the kinetic energy ek = 12 mi v 2 + μ B and toroidal canonical momentum p φ = qi ψ − mi I v  / B of the particles should be conserved during the

orbit integration. Two test particle orbits are integrated using the 2nd order Runge–Kutta scheme with t = 2Ωi−1 , and the resulting orbits are plotted in Fig. B.14(a). Initially, we set v  / v = 0.7 and v  / v = 0.4 for the passing and trapped case, respectively. Both cases have the same initial kinetic energy ek = 4.5 keV. Figs. B.14(b) and B.14(c) show the time variations of the relative errors of ek and p φ for the passing and trapped case, respectively. We can see that the employed grid systems and interpolations are accurate enough to provide stable orbit integration with good conservations of ek and p φ over the long simulation time.

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

533

Fig. B.14. (a) Test particle orbits in a shaped equilibrium. Time variations of the relative errors of the kinetic energy and toroidal canonical momentum are shown in (b) and (c) for the passing and trapped case, respectively.

The fluctuating potential Φ and the spatial part of gyrokinetic distribution function f are represented in the coordinate system (χ , φ, θ) and a grid system is generated as

χmin = χ0 < χ1 < · · · < χN χ = χmax , 0 = φ0 < φ1 < · · · < φ N φ = φmax , 0 = θ0 < θ1 < · · · < θ N θ = 2π ,

χ i = χ0 + i  χ , φk = kφ,

θ j = j θ,

(B.8) (B.9) (B.10)

where χ = (χmax − χmin )/ N χ , θ = 2π / N θ , φ = φmax / N φ . The toroidal domain can be set to simulate only a part of N w toroidal wedges, for which we set φmax = 2π / N w . On this uniform grid, the derivative operators appearing in Eqs. (10) and (11) are discretized by the 4th order accurate finite difference formulas, for example in the radial direction as

Φ i −2 − 8Φ i −1 + 8Φ i +1 − Φ i +2 , 12χ −Φi −2 + 16Φi −1 − 30Φi + 16Φi +1 − Φi +2 Φi = . 12χ 2 Φi =

(B.11) (B.12)

The integral operator for the flux surface average in Eq. (10) is discretized by the Simpson rule. The numerical discretization of Eq. (10) gives a non-sparse 2D matrix, and we employ a direct matrix solver UMFPACK [28] to solve this. On the other hand, Eq. (11) results in a sparse 3D matrix as intended, and we employ an iterative matrix solver HYPRE [29] to solve the resulting matrix. For an axisymmetric equilibrium with concentric circular cross section, we can obtain an analytic solution of the gyrokinetic Poisson equation by assuming constant plasma parameters and the cylindrical limit a0 / R 0 → 0. With these assumptions, Eq. (11) can be reduced to

 2 

ρ 1 ∂ ∂ 1 1 ∂2 ˜ + β Φ˜ = s, − 2ti r Φ˜ + 2 2 Φ λDi r ∂ r ∂ r λDe r ∂θ 2

(B.13)

where s is introduced to simply represent the source terms in the right hand side of the equation. If we express the source term s using the cylindrical Bessel functions J m (r ),

s(r , θ) =

m max 

lmax 

m=−mmax l=0

sˆlm J m (αlm r )e imθ ,

(B.14)

534

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

Fig. B.15. (a) The solution for the test of the gyrokinetic Poisson solver, and (b) the decrease of the numerical errors for increasing poloidal grid number N θ = 64, 128, 256, 512. The corresponding radial grid numbers are set as N r = 24, 48, 96, 192.

where

αlm denotes the lth-zero of J m , the solution of the reduced equation is given as Φ˜ A (r , θ) =

m max 

lmax 

m=−mmax l=0

sˆlm ρti2

λ2Di

1

α + λ2 2 lm

J m (αlm r )e imθ ,

(B.15)

De

with the boundary condition ∂ Φ˜ A (0)/∂ r = 0 and Φ˜ A (1) = 0. We can use this analytic solution to test the accuracy of the numerical solver for the gyrokinetic Poisson equation, and Fig. B.15 shows the results. Here, we chose a source form with mmax = 8 and lmax = 10 and the resulting analytic solution is shown in Fig. B.15(a). With this source and solution, we performed a set of tests with different radial and poloidal grid resolutions. The error of the numerical solution Φ˜ can be estimated as



ε=

˜ rdrdθ|Φ˜ A − Φ|,

(B.16)

and Fig. B.15(b) shows the decrease of the numerical errors for increasing grid number. As expected, the numerical errors confirm that the gyrokinetic Poisson solver is 4th order accurate. The spatial structures of fluctuating eddies in strongly magnetized plasmas are highly anisotropic and elongated along the equilibrium magnetic fields. For efficient representation of the fluctuating potential and the distribution function, we employ a field-aligned interpolation technique. More specifically, for the fluctuating potential Φ , we use the procedure as follows. For a given (χ , φ, θ), we first find the radial and toroidal grid cell [χi −1 , χi ] × [φk−1 , φk ] such that χi −1 ≤ χ ≤ χi and φk−1 ≤ φ ≤ φk . Then, we calculate the radial and toroidal coefficients for the local cubic interpolation as Eqs. (B.4)–(B.7).

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

535

Next, on each flux surface χi ∈{i −2,i −1,i ,i +1} , we integrate the following equation to find the foot points of the equilibrium magnetic field line at φ = φk ∈{k−2,k−1,k,k+1} starting from (φ, θ) (see Fig. 3).



φk Θi k (θ, φ) = θ +

dφ

B · ∇θ 

B · ∇φ χ =χ

(B.17)

i

φ

For each foot point Θi k , we find the poloidal cell containing the point as θ j i k −1 ≤ Θi k ≤ θ j i k , and calculate the following interpolation coefficients.

c −2 ( y i k ) = c −1 ( y i k ) = c 0 ( y i k ) = c 1 ( y i k ) =

y i k − y 3i k

(B.18)

6 y 3i k + y 2i k − 2 y i k

2 + y i k

(B.19)

2 − 2 y 2i k − y 3i k

(B.20)

2 y 3i k + 3 y 2i k + 2 y i k

(B.21)

6

Here y i k = [Θi k (θ, φ) − θ j i k ]/θ . Then finally, we obtain the following interpolation of Φ as

Φ(χ , φ, θ) =

i +1 

k +1 

j i k +1



i =i −2 k =k−2 j = j i k −2

c i −i (x)ck −k ( z)c j − j i k ( y i k )Φi , j ,k ,

(B.22)

where x = (χ − χi )/χ and z = (φ − φk )/φ . The spatial derivatives of Φ can be calculated straightforwardly by applying the chain rule on the known functional forms, c −2 , c −1 , c 0 , c 1 , and Θi k .

 1 ∂ Φ(χ , φ, θ) = c i −i (x)ck −k ( z)c j − j i k ( y i k )Φi , j ,k ∂χ  x i , j ,k   1 ∂ ∂ 1 c i −i (x)ck −k ( z)c j − j i k ( y i k ) − c i −i (x)ck −k ( z)c j − j ( y i k ) Φ(χ , φ, θ) = Θi k Φi , j ,k i k ∂φ φ θ ∂φ i , j ,k

 1 ∂ c i −i (x)ck −k ( z)c j − j ( y i k )Φi , j ,k Φ(χ , φ, θ) = i k ∂θ θ i , j ,k

We performed sets of linear simulations with a fixed toroidal mode number (n = 14) to test the advantages of the magnetic field-aligned interpolation as compared to the conventional 3D interpolation constructed by the tensor product of the 1D interpolations in each direction. First, we compared necessary toroidal and poloidal grid resolutions to obtain the correct linear growth of the ITG instability with the toroidal mode number n = 14. Fig. B.16 shows a comparison of the fluctuating potential patterns on a flux surface obtained by two different interpolation methods. We can clearly see that the B-aligned interpolation (Fig. B.16(a)) captures the potential patterns aligned to the equilibrium magnetic field with a modest grid resolution N θ × N φ = 256 × 8. Fig. B.16(b) shows the contours of the potential values at the grid points in the conventional interpolation case with N θ × N φ = 256 × 32. In this figure, the lack of grid resolution is apparent, though the Lagrangian interpolation guarantees the continuity of interpolated potential values for the simulation. Then, these differences affect physical results of the simulation e.g. the linear growth rate. Fig. B.17 compares the growth rate of n = 14 linear ITG instability obtained with different interpolation methods and grid resolutions. If the conventional method is used, we need substantially higher grid resolutions to capture the correct linear instability, and if the grid resolution is too low, the instability does not grow as the magenta curve shows in the figure. Next, we tested the effects of the B-aligned interpolation on the time discretization and integration. We performed two sets of n = 14 linear ITG simulations with different time step sizes and interpolation methods. For the B-aligned interpolation cases, we fixed the grid resolution as N θ × N φ = 256 × 8 and varied the time step size as t ≈ 3Ωi−1 , 6Ωi−1 , . . . . For the conventional interpolation cases, the grid resolution was fixed as N θ × N φ = 512 × 32 and the time step size was

varied similarly t ≈ 3Ωi−1 , 6Ωi−1 , . . . . Fig. B.18 shows the required number of fixed-point iterations to solve the nonlinear equations for the BSL time integration with varying time step sizes. Here, we set the tolerance error level of the iterations as ε ≡ |( f n(l+) 1 − f n(l+−11) )/ f n(l+−11) | = 10−7 . If we use the B-aligned interpolation, the fixed-point iterations converge to a numerical solution even for a large time step size t ≈ 18Ωi−1 . However, in the conventional interpolation cases, the fixed-point iter-

ations fail to converge for a smaller time step size t ≈ 12Ωi−1 . In this case, even with the apparent drop of the iteration error ε < 10−6 , the growth of the potential shows an unphysical behavior (see the magenta curve in Fig. B.19).

536

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

Fig. B.16. Patterns of fluctuating electrostatic potential Φ˜ with n = 14 on a flux surface obtain by (a) the B-aligned interpolation and (b) the conventional interpolation method.

Fig. B.17. The linear growth of fluctuating electrostatic potential with a fixed toroidal mode number n = 14. Results from different grid resolutions and interpolation methods are compared.

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

537

Fig. B.18. (a) The number of fixed-point iterations required to solve the nonlinear equations for the BSL time integration, and (b) the variation of the relative error

l−1) ε ≡ |( f n(l+) 1 − f n(l+−11) )/ f n(+ 1 | as the iteration proceeds.

Generally, the required number of iterations increases as the time step size increases in both interpolation cases. However, we can clearly see that the B-aligned interpolation guarantees more stable fixed-point iterations with a wider range of convergent time step size, and therefore provides us a greater flexibility in the choice of the simulation parameters. Also, for a small enough time step size guaranteeing the convergence of both interpolation methods, the B-aligned interpolation gives a better accuracy as we can confirm by comparing the cases with t ≈ 3Ωi−1 in Fig. B.18(b). Here, the B-aligned case shows an order smaller iteration error for the same number of iterations l = 2. Fig. B.19 compares the growth rates of the linear instability obtained with different time step sizes and interpolation methods. We can see that all convergent cases with the B-aligned interpolation give the correct instability behavior, though the biggest t case requires a substantially large number of iterations l = 21. The conventional interpolation case with t ≈ 9Ωi−1 shows a slightly higher instability level γ ≈ 0.12 than the other convergent cases γ ≈ 0.11. In summary, the magnetic field-aligned interpolation can well reproduce the fluctuating components of the electrostatic potential and distribution function with modest grid resolutions. Also, it allows us to choose bigger time step sizes and provides more stable field-point iterations as compared to the conventional interpolation method. Appendix C. Conservation properties The present BSL scheme is non-conservative i.e. it does not guarantee the conservations of key physical quantities such as total mass, momentum, and energy. Their temporal changes depend on details of spatio-temporal discretization, interpolation, and simulation parameters. In the absence of source and sink term, the gyrokinetic Vlasov equation satisfies the following conservation relations [8,9]:

d dt

M tot (t ) = −Ξb [mi f ],

(C.1)

538

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

Fig. B.19. The linear growths of fluctuating electrostatic potential with n = 14. (a) In the B-aligned interpolation cases, all convergent ones give the correct instability growth rate γ ≈ 0.11. (b) However, in the conventional cases, we need smaller time step sizes (t ≈ 3Ωi−1 ∼ 6Ωi−1 ) to reproduce the correct growth rate.

d dt d dt

 L tot (t ) =



   ∂ I ˙ Φ f − Ξb mi v  f , dXdv ∇ψ · X − qi c ∂φ B qi

(C.2)

E tot (t ) = −Ξb [ H f ],

(C.3)

where the total mass (M tot ), momentum (L tot ), and energy (E tot ) of an annular volume enclosed by defined as

χ = [χmin , χmax ] are



M tot ≡

dXdvmi f ,



L tot ≡

dXdvmi



 E tot ≡

dXdv

I B

1 2

(C.4)

v f ,

(C.5)



mi v 2 + μ B

 f +

dx

1 8π



ρti2 λ2Di

|∇⊥ Φ|2 +

1

λ2De

Φ − Φ f

2

 (C.6)

,

and H ≡ 12 mi v 2 + μ B + qi Φα . Here, we introduced the symbol Ξb [G ] to represent the following boundary integral of a physical quantity G:



Ξb [G ] ≡

dθ dφ dv  dμdα



J x B ∗ (∇ χ

χ =χ1 · X˙ )G − χ =χ0





dχ dθ dφ dμdα J x B ∗ v˙  G

where J x denotes the Jacobian of the spatial coordinate (χ , φ, θ).

 v  =+ v ,max v  =− v ,max

,

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

539

Fig. C.20. Time variations of the relative errors in the total mass (red), momentum (green), and energy (blue). The errors are defined as [ M tot (t ) − M tot (0) +











˙ − qi ∂Φ ) f ]/ M tot (0) R 0 v ti , [ E tot (t ) − E tot (0) + dt Ξb ( H f )]/ E tot (0), respectively. dt Ξb (mi f )]/ M tot (0), [ L tot (t ) − L tot (0) + dt Ξb ( i B  f ) − dt dXdv( ci ∇ψ · X ∂φ (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) m Iv

q

Fig. C.21. Time variations of the fluctuating potential and the radial shearing of E r at the mid-minor radius.

We performed a test nonlinear simulation to check the conservation properties of the new gyrokinetic code by turningoff all the source and sink terms in the gyrokinetic equation. With a canonical Maxwellian distribution as the initial condition in the concentric circular geometry, we set the simulation grid as N χ × N θ × N φ × N v  × N u = 192 × 256 × 32 × 128 × 8,

and t = 3Ωi−1 . The ion temperature and density profile are chosen to give R / L T i = 6.92 and ηi = 3.1 at the mid-minor radius. Other geometric parameters are chosen as the same with the circular case in Section 4. The simulation domain is set as [χmin , χmax ] = [0.01, 1.0], and the radial domain for the conservation integral is chosen as [χ0 , χ1 ] = [0.25, 0.75] to avoid artificial numerical effects from the simulation boundaries. The parallel velocity domain is set with v ,max = 5v ti . Figs. C.20 and C.21 show the simulation results. In Fig. C.20, we can see that the conservation errors increase rather slowly in the initial quasilinear phase (t < 60). Entering the nonlinear phase (t > 60), zonal flows are generated and contribute to the saturation of the turbulence as can be seen in Fig. C.21. In this nonlinear saturation process, fine scale turbulence patterns are generated by the zonal flow shearing. The overlap of this saturation period and the time interval of the rapid increases of the conservation errors (60 < t < 80) suggest that the errors are caused by the lack of grid resolution to capture the fine scale phase space structures generated by the interactions of zonal flow and turbulence. References [1] A.J. Brizard, T.S. Hahm, Rev. Mod. Phys. 79 (2007) 421. [2] W.W. Lee, J. Comput. Phys. 72 (1987) 243. [3] Z. Lin, T.S. Hahm, W.W. Lee, W.M. Tang, R.B. White, Science 281 (1998) 556.

540

[4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35]

J.-M. Kwon et al. / Journal of Computational Physics 283 (2015) 518–540

S.E. Parker, C. Kim, Y. Chen, Phys. Plasmas 6 (1999) 1709. F. Jenko, W. Dorland, M. Kotschenreuther, B.N. Rogers, Phys. Plasmas 7 (2000) 1904–1910. J. Candy, R.E. Waltz, J. Comput. Phys. 186 (2003) 545–581. X. Garbet, Y. Sarazin, V. Grandgirard, G. Dif-Pradalier, G. Darmet, Ph. Ghendrih, P. Angelino, P. Bertrand, N. Besse, E. Gravier, P. Morel, E. Sonnendrücker, N. Crouseilles, J.-M. Dischler, G. Latu, E. Violard, M. Brunetti, S. Brunner, X. Lapillonne, T.-M. Tran, L. Villard, M. Boulet, Nucl. Fusion 47 (2007) 1206. Y. Idomura, H. Urano, N. Aiba, S. Tokuda, Nucl. Fusion 49 (2009) 065029. J. Abiteboul, X. Garbet, V. Grandgirard, S.J. Allfrey, Ph. Ghendrih, G. Latu, Y. Sarazin, A. Strugarek, Phys. Plasmas 18 (2011) 082503. S. Ku, C.S. Chang, P.H. Diamond, Nucl. Fusion 49 (2009) 115021. M.J. Pueschel, T. Dannert, F. Jenko, Comput. Phys. Commun. 181 (2010) 1428–1437. S. Vadlamani, S.E. Parker, Y. Chen, C. Kim, Comput. Phys. Commun. 164 (2004) 209–213. V. Grandgirard, Y. Sarazin, X. Garbet, G. Dif-Pradalier, Ph. Ghendrih, N. Crouseilles, G. Latu, E. Sonnendrücker, N. Besse, P. Bertrand, Commun. Nonlinear Sci. Numer. Simul. 13 (2008) 81–87. N. Crouseilles, T. Respaud, E. Sonnendrücker, Comput. Phys. Commun. 180 (2009) 1730–1745. T.S. Hahm, Phys. Fluids 28 (1988) 2670. T.D. Arber, R.G.L. Vann, J. Comput. Phys. 180 (2002) 339. F. Filbet, E. Sonnendrücker, Comput. Phys. Commun. 150 (2003) 247. A. Staniforth, J. Cote, Mon. Weather Rev. 119 (1991) 2206. E. Sonnendrücker, J. Roche, P. Bertrand, A. Ghizzo, J. Comput. Phys. 149 (1999) 201–220. A. Burckel, O. Sauter, C. Angioni, J. Candy, E. Fabel, X. Lapillonne, J. Psychoanal. Cult. Soc. 260 (2010) 012006. L. Villard, P. Angelino, A. Bottino, S. Brunner, S. Jolliet, B.B. McMillan, T.M. Tran, T. Vernay, Plasma Phys. Control. Fusion 55 (2013) 074017. Y. Xiao, P.J. Catto, Phys. Plasmas 13 (2006) 082307. A. Marinoni, S. Brunner, Y. Camenen, S. Coda, J.P. Graves, X. Lapillonne, A. Pochelon, O. Sauter, L. Villard, the TCV Team, Plasma Phys. Control. Fusion 51 (2009) 055016. P.H. Diamond, C.J. McDevitt, Ö.D. Gürcan, T.S. Hahm, V. Naulin, Phys. Plasmas 15 (2008) 012303. Y. Camenen, A.G. Peeters, C. Angioni, F.J. Casson, W.A. Hornsby, A.P. Snodin, D. Strintzi, Phys. Rev. Lett. 102 (2009) 125001. https://fusion.gat.com/theory/Efitgeqdsk. N. Crouseilles, M. Mehrenberger, H. Sellama, Commun. Comput. Phys. 8 (2010) 484–510. T.A. Davies, ACM Trans. Math. Softw. 30 (2004) 165. R.D. Falgout, U.M. Yang, in: Computational Science ICCS 2002, in: Lecture Notes in Computer Science, vol. 2331, 2002, p. 632. M. Rosenbluth, F. Hinton, Phys. Rev. Lett. 80 (1998) 724. A.M. Dimits, G. Bateman, M.A. Beer, B.I. Cohen, W. Dorland, G.W. Hammett, C. Kim, J.E. Kinsey, M. Kotschenreuther, A.H. Kritz, L.L. Lao, J. Mandrekas, W.M. Nevins, S.E. Parker, A.J. Redd, D.E. Shumaker, R. Sydora, J. Weiland, Phys. Plasmas 7 (2000) 969. T.H. Stix, Phys. Fluids 16 (1973) 1260. G. Dif-Pradalier, V. Grandgirard, Y. Sarazin, X. Garbet, Ph. Ghendrih, Commun. Nonlinear Sci. Numer. Simul. 13 (2008) 65–71. Y. Idomura, S. Tokuda, Y. Kishimoto, Nucl. Fusion 43 (2003) 234. J.M. Kwon, S. Yi, T. Rhee, P.H. Diamond, K. Miki, T.S. Hahm, J.Y. Kim, O.D. Gurcan, C. McDevitt, Nucl. Fusion 52 (2012) 013004.