The δf particle-in-cell code GEM is used to study the transport “shortfall” problem of gyrokinetic simulations. In local simulations, the GEM results confirm the previously reported simulation results of DIII-D [Holland et al., Phys. Plasmas 16, 052301 (2009)] and Alcator C-Mod [Howard et al., Nucl. Fusion 53, 123011 (2013)] tokamaks with the continuum code GYRO. Namely, for DIII-D the simulations closely predict the ion heat flux at the core, while substantially underpredict transport towards the edge; while for Alcator C-Mod, the simulations show agreement with the experimental values of ion heat flux, at least within the range of experimental error. Global simulations are carried out for DIII-D L-mode plasmas to study the effect of edge turbulence on the outer core ion heat transport. The edge turbulence enhances the outer core ion heat transport through turbulence spreading. However, this edge turbulence spreading effect is not enough to explain the transport underprediction.

Gyrokinetic simulations are widely used to study tokamak plasmas and the models are now comprehensive enough that the simulation predictions can be directly compared to experimental results. Despite the success of the validation and verification studies using gyrokinetics, there are still open questions remain, such as the “transport shortfall”.1,2 In GYRO3 simulations of a DIII-D L-mode discharge dominated by only micro-instabilities, i.e., no magnetohydrodynamic (MHD) instabilities in the model, the ion heat flux transport, Qi, from simulation agrees well with the experimental result predicted by the transport code ONETWO at the inner core for ρ=0.5, where ρ is the square root of the normalized toroidal flux. However, the simulation Qi is much smaller than the experimental result at outer core, for ρ=0.75. It seems this underprediction cannot be explained by the uncertainties of the experimental parameters used in the simulations. The problem has even raised questions on the validity of the gyrokinetic model,4,5 and also raised questions regarding whether ONETWO is a good model for experimental transport predictions.6 

On the other hand, a GYRO simulation of the outer core region of an Alcator C-Mod L-mode observes underprediction for Qe but not for Qi.7 In another systematic transport analysis of Alcator C-Mod L-mode using GYRO at different radial locations for two discharges differing in terms of the amount of auxiliary ion heating (termed as low and high power discharges),8,9 no DIII-D L-mode type Qi shortfall is found whatsoever and Qe underprediction is observed in the case of the low power discharge. This Qe underprediction is likely due to the lack of the high k electron scale physics in the simulations. A more recent multi-scale simulation including electron temperature gradient (ETG) scale fluctuations for these cases10,11 finds the Qe underprediction is diminished. Furthermore, simulations of the L-mode plasmas in the ASDEX Upgrade tokamak using the nonlinear gyrokinetic code GENE12,13 have shown good agreement with experimental results and there is no shortfall within the uncertainties of experimental profiles in both the inner and outer core regions.14 Then, a new study of the same DIII-D shortfall case in Refs. 1 and 2 using the GENE code shows no shortfall by varying the pressure profiles.15 

In this paper, we address the shortfall problem with another gyrokinetic code, GEM.16,17 Unlike the GYRO and GENE code which are both continuum (Eulerian) models, GEM is a δf particle-in-cell code, and it has been well benchmarked with continuum codes such as GYRO and GS2 in different studies.18,19 A particle code provides a critical verification of the previously reported nonlinear calculations using a completely different numerical approach. In Ref. 2, an earlier version of GEM has been used to simulate the same DIII-D shortfall case and has similar results as GYRO. An alternative flux-tube model has been implemented recently using a complete local equilibrium solution of the Grad-Shafranov equation,19 and we will apply this model in the current study.

Additionally, the GYRO and GENE models used in previous shortfall studies are local, i.e., flux tube simulations and the heat flux transport might be affected by this restriction. Adding edge physics, such as turbulence spreading20–24 and avalanches propagating from the edge,25 could increase the transport in the outer core region. GEM is an intrinsically global code and it has been successfully used to study both core26–28 and edge physics.29–31 Here, in addition to flux tube simulations, we also apply the global GEM model to study the edge turbulence spreading effects on the outer core transport.

The purpose of this paper is to document a fairly comprehensive GEM study of the relevant shortfall cases, and explore the mechanism of global edge turbulence spreading that has the potential to explain the shortfall. We first analyze the L-mode discharges from DIII-D and Alcator C-Mod tokamaks of previous GYRO studies in Refs. 1, 2, 8, 9, and 32, and compare ion and electron heat fluxes obtained from our simulations with those estimated in experimental measurements. In addition to these nonlinear results, we will identify the dominant linear instabilities at different radial locations across the minor radius of the two tokamaks. A global calculation is then carried out to examine the global profile effects and the effect of edge turbulence on the shortfall for the DIII-D L-mode discharge. We find that GEM has reconfirmed the GYRO results, i.e., the shortfall exists in the DIII-D case and there is no DIII-D L-mode like shortfall in Alcator C-Mod case, and although the edge physics has significant effects on the outer core transport, the effect is not big enough to explain the shortfall.

In this section, we document the δf particle simulation model in the GEM code. The main features include electromagnetics with full kinetic electron physics, δf particle in cell method, p formulation, a split weight scheme for electrons, high β Ampere algorithm, a field line following coordinate system covering 0θ<2π and global profile effects. The model includes collisions, equilibrium flow, an arbitrary shaped tokamak equilibrium, and impurity species. In terms of the canonical momentum ps=vs+qsmsA, where s stands for particle species, the gyrokinetic equation is given as

fst+vGs·fs+ṗsfsp=C(fs),
(1)

where,

ṗ,s=qsmsb̃·ϕμsmsb̃·B+v,s(b·b)·vE+qsmsvGs·A.

The guiding center velocity is given as vGs=vsb̃+vds+vE, where b̃=b+δB/B,vE=E×b/B, and for low β plasmas

vds=msv2/2qsB3B×B+msv2qsBb×(b·b),=ms(v2+v2/2)qsB3B×B+msv2qsB3b×[(×B)×B],ms(v2+v2/2)qsB3B×B.
(2)

Note that in Eq. (1), C(fs) is the collision operator. It is considered that ion collision operator C(fi)=0. The Lorentzian operator is considered for electrons so that C(fe)=CL(fe), where

CL(fe)=0.5νeλ(1λ2)λfe,

with λ=v/v is the pitch angle. Also

νe=n0ee4lnΛ4πϵ02me2v3(Zeff+Heemev2/2T0e),

where

Hee(x)=ex2πx+(112x2)erf(x).

The equilibrium distribution function is considered to be Maxwellian for all species, given by, f0s(x)=n0s(x)/((2π)1.5vTs3(x))eεs/msvTs2, where εs=ms(vs2+ps2)/2 is the total energy for species, s. Note that εs is a function of ps instead of vs.

A δf scheme is considered for ions, for which the evolution can be expressed as

dδfidt=(viδBB+vE)·f0iε̇if0iεi.

A split weight scheme33,34 is implemented for electrons, where the adiabatic component of the perturbed distribution function is split into two parts, as fe=f0eϵgϕ(f0e/εe)+h. Then the perturbed electron distribution function evolves as

dhdtCL(fe)=(veδBB+vE)·f0eε̇ef0eεe+ϵg(ϕt+vGe·ϕ)f0eεe.
(3)

In the last two equations ε̇s=μsvGs·B+mspsṗs. The weight equations for ions and electrons wi and we can be derived from the evolution equations for δfi and h and are as follows:

dwidt=(viδBB+vE)·f0if0iε̇i1f0if0iεi,
(4)
dwedt=(veδBB+vE)·f0ef0eε̇e1f0ef0eεe+ϵg(ϕt+vGe·ϕ)1f0ef0eεe+ϵgϕ[(veδBB+vE)·f0eεe+ε̇ef0eεe].
(5)

The electrostatic potential ϕ is calculated using the gyrokinetic Poisson equation

(ϕϕ̃)+ϵgτϕ=δfiδ(R+ρx)dRdvhdv,
(6)

where ϕ̃=kΓ0(k2ρti2)ϕkeik·x,ϕ=kϕkeik·x, and dv=vdvdvdξ, with ξ and ρ, respectively, are the gyroangle and the vector from particles' gyrocenter to the actual positions of the particles. The vector potential, similarly, is given by Ampere's law

(2+βeme)A=βi(δfipδ(R+ρx)dRdv(ϵgϕf0eεe+h)vdv).
(7)

It is to be noted that ϕ/t required in the calculation of dh/dt and dwe/dt in Eqs. (3) and (5) can be obtained by taking time derivative of the Poisson equation expressed in Eq. (6).

The field line following coordinate system has been used in the model. They are given as, x=rr0,y=r0q0(0θdθq̂(r,θ)ζ),z=q0R0θ, where (r,θ,ζ) are the toroidal coordinates, R0 and r0, respectively, are the major radius at the magnetic axis and the minor radius at the center of the simulations box, q0 is the safety factor at r0, and q̂(r,θ)=B·ζ/B·θ. Thus, (x,y) labels the field line while z is the coordinate along the field line. The particle motion is given by ẋ=vG·x,ẏ=vG·y, and ż=vG·z. Defining mp = proton mass and T0i = ion temperature, vu=T0i/mp,xu=vu/ωci,ωci=eB0/mp,tu=1/ωci the normalization scheme can be seen in Tables I and II. For further details of the flux tube and global model of GEM readers are referred to Refs. 16 and 17.

TABLE I.

Normalization I.

Parameter f s 0 δ f s ϕ A j
Normalization  n 0 v u 3   n 0 v u 3   T 0 i / e   T 0 i / e v u   e n 0 v u  
Parameter f s 0 δ f s ϕ A j
Normalization  n 0 v u 3   n 0 v u 3   T 0 i / e   T 0 i / e v u   e n 0 v u  
TABLE II.

Normalization II.

Parameter t x y z v ms Ts
Normalization  tu   xu   xu   xu   vu   mp   T 0 i  
Parameter t x y z v ms Ts
Normalization  tu   xu   xu   xu   vu   mp   T 0 i  

We first study the “standard” DIII-D L-mode shortfall case1,2 discharge 128913. The parameters corresponding to this case are obtained from Refs. 1 and 35 and delineated in Table III for two radial locations ρ=0.5 and ρ=0.75, in which R is the major radius and a is the minor radius, R/a is aspect ratio, κ is elongation, sκ=rdlnκ/dr, δ is triangularity, sδ=rdδ/dr, Δ is Shafranov shift, q is the safety factor at the given values of ρ, ŝ is the shear, Te and Ti are the temperatures of electron and ion, ρ=ρs/a,ρs=cs/ωi,ωi=eB0/mi, νe is the collision frequency, and βe=2μ0neTe/B2 is the electron β. The pressure gradient terms are denoted by R/Lns=Rdlnn0s/dx and R/LTs=RdlnT0s/dx for species s=e,i. γE×B is the equilibrium radial electric field shearing rate. The values of βe and νe are estimated using the relations βe=neTe/(B2/2μ0) and νe=nee4lnΛ/4πϵ02m2vthe3.

TABLE III.

Parameters for the DIII-D L-mode plasmas.

Parameters ρ = 0.5 ρ = 0.75
ne (1019m3)  2.08  1.65 
Te (keV)  0.964  0.433 
Ti (keV)  0.805  0.511 
R/a  2.81  2.79 
κ   1.30  1.36 
s ̂ κ   0.0501  0.245 
δ0   0.153  0.237 
s δ 0   0.180  0.438 
Δ  −0.0863  −0.105 
q   1.83  2.77 
s ̂ = r d ln q / d r   0.625  2.06 
T i / T e   0.835  1.18 
R / L T i   5.0861  7.1145 
R / L n e   3.0067  3.0132 
R / L T e   7.4184  13.7547 
γ E × B ( c s / a )   0.0402  0.0782 
Z eff   1.32  1.33 
ρ   0.0036  0.0024 
ν e ( c s / a )   0.112  0.44 
βe   0.0018  0.00065 
Parameters ρ = 0.5 ρ = 0.75
ne (1019m3)  2.08  1.65 
Te (keV)  0.964  0.433 
Ti (keV)  0.805  0.511 
R/a  2.81  2.79 
κ   1.30  1.36 
s ̂ κ   0.0501  0.245 
δ0   0.153  0.237 
s δ 0   0.180  0.438 
Δ  −0.0863  −0.105 
q   1.83  2.77 
s ̂ = r d ln q / d r   0.625  2.06 
T i / T e   0.835  1.18 
R / L T i   5.0861  7.1145 
R / L n e   3.0067  3.0132 
R / L T e   7.4184  13.7547 
γ E × B ( c s / a )   0.0402  0.0782 
Z eff   1.32  1.33 
ρ   0.0036  0.0024 
ν e ( c s / a )   0.112  0.44 
βe   0.0018  0.00065 

Despite the extensive nonlinear analysis of these two radial locations in previous studies that have found the ion heat flux transport shortfall at ρ=0.75 but no shortfall at ρ=0.5, a complete linear analysis is not yet available. Here, we first run linear flux tube simulations to identify the dominant instabilities at the two locations. For the sake of simplicity, we remove the effect of Er shearing rate in these linear simulations. The various linear scans are discussed below in tandem.

Figure 1 shows the dispersion relations obtained from linear GEM simulations for ρ=0.5 and ρ=0.75. The real frequency and linear growth rate of the dominant instability are plotted as functions of kθρs. The two locations have similar results, and we identify the dominant instability to be the ion temperature gradient (ITG) mode. The real frequencies in Figs. 1(a) and 1(c) are negative, corresponding to the ion diamagnetic direction in GEM's convention. The corresponding growth rates are shown in Figs. 1(b) and 1(d). Furthermore, a βe scan shows that the mode is electrostatic, and a collision frequency scan shows the mode is independent of collision frequency. These are all typical ITG characteristics.

FIG. 1.

Linear dominant instability for DIII-D L-mode case: (a) real frequency, (b) growth rate for ρ=0.5; (c) real frequency, (d) growth rate for ρ=0.75 plotted varying kθρs.

FIG. 1.

Linear dominant instability for DIII-D L-mode case: (a) real frequency, (b) growth rate for ρ=0.5; (c) real frequency, (d) growth rate for ρ=0.75 plotted varying kθρs.

Close modal

For a further confirmation, we can identify the main drive that provides free energy to these modes at the two locations. We carry out electron and ion temperature gradient scans as represented by R/LTe and R/LTi. In Fig. 2, the electron temperature gradient is varied by 10% around the experimental value and apparently the growth rates are insensitive to R/LTe at both locations. In the ion temperature gradient scan shown in Fig. 3, the growth rates increase linearly with R/LTi. Therefore, the mode is indeed driven by the ion temperature gradient.

FIG. 2.

R/LTe scan for growth rate at ρ=0.5, and 0.75 for kθρs=0.7 and kθρs=0.66, respectively.

FIG. 2.

R/LTe scan for growth rate at ρ=0.5, and 0.75 for kθρs=0.7 and kθρs=0.66, respectively.

Close modal
FIG. 3.

R/LTi scan for growth rate at ρ=0.5, and 0.75 for kθρs=0.7 and kθρs=0.66, respectively.

FIG. 3.

R/LTi scan for growth rate at ρ=0.5, and 0.75 for kθρs=0.7 and kθρs=0.66, respectively.

Close modal

Nonlinear flux tube simulation results of the ion heat flux transport of the two locations are shown in Fig. 4. The Qi values are time-averaged over the steady state after saturation that included time window of width at least 1.5×105 in tωci. The experimental values (circles) are obtained from Ref. 2 for comparison. The simulation results using the “standard” experimental parameters are smaller than the ONETWO prediction of the experimental results at both locations. At ρ=0.5, the simulation result is about half of the experimental value, while at ρ=0.75, the simulation result is much smaller, and hence, the shortfall. We note here that our simulation result at ρ=0.5 is smaller than the old GEM results presented in Ref. 2, this is because we have updated the flux tube simulation algorithm and collision operator since then. The general trend that ρ=0.75 has a significant shortfall, however, is consistent with previous studies.1,2

FIG. 4.

Nonlinear simulation results of Qi with the original parameters (star) and varied R/LTi values compared to the experimental results (circle) at ρ=0.5 and ρ=0.75.

FIG. 4.

Nonlinear simulation results of Qi with the original parameters (star) and varied R/LTi values compared to the experimental results (circle) at ρ=0.5 and ρ=0.75.

Close modal

Could this under-prediction by simulation be due to uncertainties of the measured plasma parameters? Since the dominant instability is ITG, the Qi flux is most sensitive to ion temperature gradient. Here, we increase the value of R/LTi by 10%, 25%, 30%, and 50% to test the dependence of ion heat transport on R/LTi. These new results are also shown in the same figure. At ρ=0.5, the simulation value matches the experimental one when R/LTi is increased by 25% (diamond); while at ρ=0.75, the simulation result for ion heat transport agrees with the experimental value when R/LTi is increased by 30% (square). These ranges are apparently beyond the error of the experimental measurement.

Next, we consider the two Alcator C-Mod cases studied in Refs. 8 and 9. The two cases are labeled as “low power” (LP) and “high power” (HP) characterized by the amount of auxiliary heating of ions. The physical parameters obtained from Ref. 8 are listed in Tables IV and V for two locations, r/a=0.5 and r/a=0.8, where r is the Miller parameter that labels the flux surface and a is the minor radius. We note that in Ref. 8 the driving gradients and suppression terms are slightly adjusted to match the heat flux from GYRO simulation with the experimental heat flux obtained from TRANSP36 within the experimental uncertainty, and thus the simulation study is termed as the “flux matched simulation.”8 The GEM flux tube simulations here use these same parameters.

TABLE IV.

Parameters for Alcator C-Mod case.

Power r/a R/a κ s κ δ s δ Δ q0 s ̂ Te/Ti TB/Ti Z eff ρ ν e ( c s / a ) βe
LP  0.5  3.10  1.23  0.06  0.09  0.11  −0.05  1.30  0.91  1.124  1.0  1.90  0.0036  0.11  2.24e−3 
LP  0.8  3.08  1.33  0.36  0.21  0.55  −0.09  2.64  2.33  0.855  1.0  1.90  0.0017  0.52  5.66e-4 
HP  0.5  3.12  1.23  0.07  0.09  0.12  −0.09  1.35  0.91  1.37  1.0  2.75  0.0049  0.03  4.21e-3 
HP  0.8  3.09  1.33  0.36  0.22  0.57  −0.15  2.70  2.34  0.935  1.0  2.75  0.0025  0.14  1.2e-3 
Power r/a R/a κ s κ δ s δ Δ q0 s ̂ Te/Ti TB/Ti Z eff ρ ν e ( c s / a ) βe
LP  0.5  3.10  1.23  0.06  0.09  0.11  −0.05  1.30  0.91  1.124  1.0  1.90  0.0036  0.11  2.24e−3 
LP  0.8  3.08  1.33  0.36  0.21  0.55  −0.09  2.64  2.33  0.855  1.0  1.90  0.0017  0.52  5.66e-4 
HP  0.5  3.12  1.23  0.07  0.09  0.12  −0.09  1.35  0.91  1.37  1.0  2.75  0.0049  0.03  4.21e-3 
HP  0.8  3.09  1.33  0.36  0.22  0.57  −0.15  2.70  2.34  0.935  1.0  2.75  0.0025  0.14  1.2e-3 
TABLE V.

Parameters for Alcator C-Mod case.

power r/a n e ( 10 20 ) Te (keV) R/Lni R/LTi R/Lne R/LTe R/LnB R/LTB n D / n e n B / n e γ E × B ( c s / a )
LP  0.5  1.15  1.41  2.356  5.1460  2.356  8.339  2.356  5.1460  0.875  0.025  0.005 
LP  0.8  0.79  0.52  6.006  11.2112  6.006  13.3672  6.006  11.2112  0.875  0.025  −0.023 
HP  0.5  1.18  2.59  2.7144  3.9000  2.7144  7.6128  2.7144  3.9000  0.90  0.02  −0.013 
HP  0.8  0.82  1.05  4.8513  6.7053  4.8513  9.9189  4.8513  6.7053  0.90  0.02  0.004 
power r/a n e ( 10 20 ) Te (keV) R/Lni R/LTi R/Lne R/LTe R/LnB R/LTB n D / n e n B / n e γ E × B ( c s / a )
LP  0.5  1.15  1.41  2.356  5.1460  2.356  8.339  2.356  5.1460  0.875  0.025  0.005 
LP  0.8  0.79  0.52  6.006  11.2112  6.006  13.3672  6.006  11.2112  0.875  0.025  −0.023 
HP  0.5  1.18  2.59  2.7144  3.9000  2.7144  7.6128  2.7144  3.9000  0.90  0.02  −0.013 
HP  0.8  0.82  1.05  4.8513  6.7053  4.8513  9.9189  4.8513  6.7053  0.90  0.02  0.004 

For flux tube linear analysis, we first carry out a kθρs scan for real frequency and growth rate at the two locations for both the low power and high power cases without Er. It appears the low power case is dominated by the ITG at both locations, as shown in Fig. 5. At r/a=0.5 in Figs. 5(a) and 5(b), the real frequency is in the ion diamagnetic direction, and the growth rate increases with kθρs because of magnetic drift resonance which depends on k, then reaches a maximum value and eventually rolls off due to the finite Larmor radius effect. At r/a=0.8 as in Figs. 5(c) and 5(d), the real frequency is in the electron direction at lower k and goes into the ion direction at higher k. The growth rate continues to increase with increasing kθρs until kθρs0.7 for a really high n = 120 and tends to saturate. The ion mode is clearly dominant.

FIG. 5.

Real frequency (left) and growth rate (right) for the low power discharge case. Upper panels are for r/a=0.5, while lower panels are for r/a=0.8.

FIG. 5.

Real frequency (left) and growth rate (right) for the low power discharge case. Upper panels are for r/a=0.5, while lower panels are for r/a=0.8.

Close modal

The high power case shown in Fig. 6 has a different story. At r/a=0.5, the real frequency is first in the ion direction with a small amplitude when kθρs is smaller than 1.1, then jumps to the electron direction with a larger amplitude for bigger k and then increases with k, shown in Fig. 6(a). It appears the dominant instability here is an ion mode for low k and an electron mode for high k, and the ion mode has higher growth rate, as shown in Fig. 6(b). Since for the high power case R/LTi and νe are decreased and Te is increased, we identify this electron instability as the trapped electron mode (TEM). At r/a=0.8 shown in Figs. 6(c) and 6(d), the real frequency is in the electron diamagnetic direction at low k then switches to the ion diamagnetic direction and increase with k. The growth rates increase with kθρs, asymptoting at higher k. The absence of a strong finite Larmor radius effect could be due to the presence of trapped electron modes in this discharge. Thus, the high power discharge shows a dominant TEM not found in the low power cases. Also, radially outer locations are more susceptible to TEM activity, partly because of the increased fraction of trapped electrons as one moves from the core towards the edge of the tokamak.

FIG. 6.

Real frequency (left) and growth rate (right) for the high power discharge. Upper panels are for r/a=0.5, while lower panels are for r/a=0.8.

FIG. 6.

Real frequency (left) and growth rate (right) for the high power discharge. Upper panels are for r/a=0.5, while lower panels are for r/a=0.8.

Close modal

βe scans show the growth rate is independent of βe in all cases, indicating the modes are electrostatic. This is expected because the βe is usually very small in Alcator C-Mod plasmas. A collision scan is also carried out for the two cases. The experimental value of νe is increased and decreased by 10%. It is observed that the variation in the growth rate with collision frequency is weak around the experimental value.

To identify the drive of the instabilities, we look at the sensitivity of the mode growth rates with respect to the temperature gradients of electrons and ions. First, we vary the electron temperature gradient by ±10% keeping the value of ion temperature gradient R/LTi fixed. The corresponding growth rates are calculated and displayed in Fig. 7 for all the cases. The mode growth rate is virtually independent of R/LTe except at r/a=0.8 for the high power case, in which γ increases linearly with R/LTe, confirming that the instability here is the TEM. In the corresponding R/LTi scan as in Fig. 8, except for r/a=0.8 in the high power case which is also the TEM dominant, the dominant instabilities are driven by the ion temperature gradient for the chosen value of kθρs.

FIG. 7.

Dependence of mode growth rates on the electron temperature gradient. Upper panels are for the low power discharge at r/a=0.5 (left) and r/a=0.8 (right) for kθρs=0.66 and kθρs=0.4, respectively; while lower panels are for the high power discharge at r/a=0.5 (left) and r/a=0.8 (right) for kθρs=0.9 and kθρs=0.6, respectively.

FIG. 7.

Dependence of mode growth rates on the electron temperature gradient. Upper panels are for the low power discharge at r/a=0.5 (left) and r/a=0.8 (right) for kθρs=0.66 and kθρs=0.4, respectively; while lower panels are for the high power discharge at r/a=0.5 (left) and r/a=0.8 (right) for kθρs=0.9 and kθρs=0.6, respectively.

Close modal
FIG. 8.

Dependence of mode growth rates on ion temperature gradient. Upper panels are for the low power discharge at r/a=0.5 (left) and r/a=0.8 (right) for kθρs=0.66 and kθρs=0.4, respectively; while lower panels are for the high power case at r/a=0.5 (left) and r/a=0.8 (right) for kθρs=0.9 and kθρs=0.6, respectively.

FIG. 8.

Dependence of mode growth rates on ion temperature gradient. Upper panels are for the low power discharge at r/a=0.5 (left) and r/a=0.8 (right) for kθρs=0.66 and kθρs=0.4, respectively; while lower panels are for the high power case at r/a=0.5 (left) and r/a=0.8 (right) for kθρs=0.9 and kθρs=0.6, respectively.

Close modal

For nonlinear flux tube results, the time averaged ion and electron heat fluxes Qi and Qe from the GEM simulations (stars) are compared to the experimental values (squares) obtained from Ref. 8, shown in Fig. 9. The simulated heat fluxes via ion and electron channels are in good agreement with those from experiments evaluated using the TRANSP36 code except in the low power case, where electron heat flux is lower than the experimental range at both radial locations. This is in good agreement with GYRO results in Ref. 8.

FIG. 9.

Comparison of GEM results with experimental values for the locations r/a=0.5 and r/a=0.8. Upper panels: ion heat flux (left) and electron heat flux (right) for the low power discharge case. Lower panels: ion heat flux (left) and electron heat flux (right) for the high power discharge case. The squares show the experimental range for a given radial location, while the stars show the estimation from GEM simulations.

FIG. 9.

Comparison of GEM results with experimental values for the locations r/a=0.5 and r/a=0.8. Upper panels: ion heat flux (left) and electron heat flux (right) for the low power discharge case. Lower panels: ion heat flux (left) and electron heat flux (right) for the high power discharge case. The squares show the experimental range for a given radial location, while the stars show the estimation from GEM simulations.

Close modal

Thus, in contrast to the results from DIII-D L-mode plasmas where the ion heat flux is low (i.e., the shortfall) in the outer radial locations, the present analysis indicates no such Qi shortfall exists in these Alcator C-Mod L-mode discharges. However, the electron heat flux is under-predicted at all radial locations in the low power discharge. This could be explained by the fact that the low power discharge is characterized by higher a/LTi, higher collisions and lower electron temperatures8 and hence electron modes are less dominant, resulting in a lower electron heat flux. In fact, including high k, ρe scale resolution in simulations of the low power discharge showed better agreement between experimental electron heat flux and the simulation results.10,11 The including of high k turbulence may, in principle, also enhance the ion heat fluxes in those simulations, but the effect is within the experimental error bars.

In this section, we study the global effects on the local transport level of DIII-D discharge 128913. Flux tube simulations, though with a finite radial domain, are local studies that assume the pressure and magnetic equilibrium parameters and their gradients are homogeneous across the simulation domain. The only exception is the safety factor q, which is implemented using a linear profile matching the local values of q0 and ŝ, to maintain the magnetic shear. In GEM's global model,29 the realistic pressure and magnetic equilibrium profiles are implemented including the experimental q profile. With the global model different physical instabilities may appear at different radial locations and couple with each other within the simulation domain, and turbulence could grow and spread across the domain.

For the present study, the experimental magnetic equilibrium profile (i.e., the eqdsk file) and pressure profile are used directly in GEM, and all radial pressure gradients are then calculated numerically within the simulation. As a result, the Miller parameterization37 and other simulation parameters used in our global simulations are slightly different from the “standard parameters” used in other studies of the shortfall problem.1,2 We first compare flux tube GEM simulation results using these parameters to the experimental prediction of ONETWO2 in Fig. 10. Here, the simulation value of Qi at ρ=0.5 is quite close to the one we get with the standard parameters, though now smaller than the experimental prediction. At ρ=0.75 there is apparently a Qi shortfall. In addition to ρ=0.5, the simulations show a second peak of Qi at ρ=0.8, and the dominant instability is ITG at both locations. The TEM appears to dominate around ρ=0.9, with a very low Qi. At ρ=0.6 to 0.65, however, the simulations are almost completely stable.

FIG. 10.

GEM flux tube simulation results of the ion heat flux transport Qi compared to the ONETWO prediction of the experimental results.

FIG. 10.

GEM flux tube simulation results of the ion heat flux transport Qi compared to the ONETWO prediction of the experimental results.

Close modal

Although the simulation results at all locations are smaller than the ONETWO prediction, the trend of instabilities is consistent with the pressure profile. Figure 11 shows the variation of ηi and ηe. ηi peaks around ρ=0.5 and ρ=0.8, the same locations of ITG dominance, and ηi is below 2 around ρ=0.65, where ITG is stable. For ρ>0.9, ηe is big but ηi is very small, and hence we have a small Qi. In contrast to the smooth experimental fluxes, the local simulation fluxes have ups-and-downs. One possible mechanism that may reconcile the local simulation and experimental results is that turbulence spreading from the strongly unstable regions could make up for the shortfall in the stable regions. Especially, strong turbulence from the steep gradient region of the edge, which is not included in any local studies, may help increase the overall transport of the outer core region. We need to test this possible mechanism with global simulations.

FIG. 11.

The ηi and ηe variation of the pressure profile used in the global simulation study.

FIG. 11.

The ηi and ηe variation of the pressure profile used in the global simulation study.

Close modal

The GEM global simulations here focus on the outer core to edge region, and cover the radial domain of 0.699r/a0.999, or 0.641ρ0.998 in terms of ρ. The simulations are electromagnetic and collisional, with gyrokinetic deuterium and carbon ions and drift kinetic electrons. The particle-in-cell scheme has 128, 128, and 32 grids in Lx, Ly, and Lz in terms of the field-line-following coordinates17 and 32 particles per species per cell. The simulation time step is Δt=1/Ωi, where Ωi is the proton gyro-frequency at the core.

Global linear results are consistent with the flux tube simulations discussed in Sec. III. The dominant instability appears at the edge steep density gradient region and is identified as the TEM. In a simulation with a smaller radial domain that includes the edge region only, we have the same linear growth rate as in the simulation with the full domain. In a simulation that does not include the edge, we have observed the ITG instability that appears in the outer core region. Also same as in the flux tube runs, high-n modes (n > 50) dominate both the ITG and the TEM. In nonlinear global simulations including the low-n modes, we have observed a geodesic acoustic mode (GAM)-like instability that has a saturation level beyond the regime where the δf method is valid, hence here we keep the n20 modes only to exclude any low-n instabilities, while the kr>5·2π/Lx zonal flow is included.

To test the edge effects on the core transport, we have two nonlinear simulations: one with the original experimental profile and the other with flattened profiles at the edge. Figure 12 shows the ne, Ti, and Te profiles for the two simulations. In the flattened case, the pressure gradient terms are multiplied by a fattening function exp[((r/a0.85)/0.05)2] for r/a>0.85, and therefore the ITG dominant region, which is inside r/a<0.85, is not affected by the flattening. Without the pressure gradient drive, the flattened case is almost linearly stable in the edge region.

FIG. 12.

The global simulation profiles for electron density, deuterium and carbon ion temperature, and electron temperature. Solid line is with the original edge, and dashed line is with the flattened edge.

FIG. 12.

The global simulation profiles for electron density, deuterium and carbon ion temperature, and electron temperature. Solid line is with the original edge, and dashed line is with the flattened edge.

Close modal

The nonlinear evolution of the two simulations is shown in Fig. 13. The root mean square (RMS) value of the overall volume-averaged perturbed electrostatic potential is plotted to indicate the turbulence level. In the case with the original edge, the linear growth is very strong and the simulation saturates very quickly at a high level. In the case with the flattened edge, the linear stage lasts for a very long time and the simulation eventually saturates at a low level. Apparently, the edge instability has significantly increased the total turbulence.

FIG. 13.

The nonlinear evolution of the volume-averaged RMS δϕ for the two global simulations.

FIG. 13.

The nonlinear evolution of the volume-averaged RMS δϕ for the two global simulations.

Close modal

Figures 14 and 15 show the process of turbulence spreading of the two cases. Snapshots of the perturbed electrostatic field δϕ(x,y) turbulence were taken from three times during the evolution. In Fig. 14 with the full edge pressure gradient drive, the TEM instability at edge dominates the whole simulation domain for a long time, as shown in the top panel for tΩi=2000. The ITG instability at 0.7<ρ<0.8 grows very slowly and then turbulence appear to peak at two locations, one at the edge from the TEM and other at the outer core from the ITG, as shown in the second panel. During this period, the edge instability has already saturated and its turbulence is spreading inward towards the core. This phenomenon is quite similar to the turbulence spreading studied in Ref. 22 in which ITG turbulence spreads from the outer core region (around r/a = 0.7) to the core. The TEM and the ITG start independently and grow locally for a while. As the ITG instability saturates, its turbulence also begins to spread both inward and outward. Turbulence from the two instabilities goes across the linearly stable region around ρ=0.85 and eventually spreads all over the simulation domain, with visible zonal flow structures, shown in the bottom panel.

FIG. 14.

Evolution of the δϕ(x,y) turbulence for the original edge case.

FIG. 14.

Evolution of the δϕ(x,y) turbulence for the original edge case.

Close modal
FIG. 15.

Evolution of the δϕ(x,y) turbulence for the flattened edge case.

FIG. 15.

Evolution of the δϕ(x,y) turbulence for the flattened edge case.

Close modal

In the case with the flattened edge pressure profile as shown in Fig. 15, the edge region is stable and the only instability is the ITG at the outer core. After a very long linear growth, the ITG turbulence spreads across the linearly stable region into the edge. Similar to Ref. 22, large amplitude turbulence is still largely concentrated in the linearly unstable region even after the spreading.

Figure 16 shows the time-averaged turbulence intensity, which is represented by RMS δϕ(x) averaged over y and z, for the two cases at the quasi steady state after saturation. In the case with the original edge pressure profile (solid line), turbulence slightly peaks at the TEM dominant edge region, but it has spread throughout the simulation domain. In the case with the outer core ITG only (dashed line), turbulence reaches the edge region but it still largely concentrates in the outer core. For everywhere including the ITG dominant region of 0.7<ρ<0.8, the turbulence intensity with the edge instability is clearly higher than that without the edge instability. Therefore, we conclude that the edge instability has increased the fluctuation level at the outer core region. This is indeed a significant global effect that it has even greatly increased the turbulence level at inner boundary of the simulation domain near ρ=0.65, which is linearly stable in local simulations.

FIG. 16.

Time averaged turbulence intensity after saturation for the two cases.

FIG. 16.

Time averaged turbulence intensity after saturation for the two cases.

Close modal

How does the turbulence spreading from the edge affect the outer core transport? In Fig. 17, we compare Qi at the saturation stage from our global simulations to the flux tube simulations and the experimental results. Results from the case of the global simulation without the edge drive qualitatively agree with that of the flux tube simulations. Compared to the flux tube results, the global results are higher for ρ<0.8 and lower for ρ>0.8, but the trend is the same. The transport peaks around ρ=0.8 and diminishes in both ends. When the edge drive is included as in the results shown in the solid line, the transport at ρ=0.75 is strongly increased to more than twice of the results from flux tube simulations, apparently an effect of turbulence spreading from edge. The TEM edge instability also increases Qi in the edge region, where the ITG is stable. The edge effect, however, is not strong enough to make up for the Qi shortfall, especially in the ITG stable regions. The carbon impurity heat flux is included in the calculation of Qi here but its contribution is very small.

FIG. 17.

Qi at the saturation stage from global simulations with the original edge pressure profile (solid line), the flattened edge profile (dashed line), flux tube simulations (squares) compared to the ONETWO predictions.

FIG. 17.

Qi at the saturation stage from global simulations with the original edge pressure profile (solid line), the flattened edge profile (dashed line), flux tube simulations (squares) compared to the ONETWO predictions.

Close modal

In Fig. 18 for Qe, edge turbulence also increases the transport at the outer core, but the increase is not as big as for Qi. With the pressure gradient drive there is a second peak of Qe at the edge, where Qe is almost zero in the case with the flattened edge.

FIG. 18.

Qe of the two global cases compared to the ONETWO predictions.

FIG. 18.

Qe of the two global cases compared to the ONETWO predictions.

Close modal

Figure 19 shows the ion profile relaxation as represented by δni and δpi. In δf simulations for f=f0+δf, we assume δff, and Fig. 19(a) indicates that this conditions is well satisfied even after a very long time on nonlinear run. Since δni/ni0 is within 1% and δpi/pi0<4%, our δf simulations are valid and this scale of relaxation is physical. Figs. 19(b) and 19(c) show that the relaxation tend to decrease the pressure gradients at the location with strongest linear instability (around ρ=0.75), and as a result increase the gradients at nearby locations. Lower pressure gradients correspond to lower instability, therefore profile relaxation is a mechanism that lowers the transport at the most unstable region. This phenomenon is consistent with the physical intuition that transport tends to smooth the particle and energy distributions. Taking into account for this effect, one would expect simulations without profile relaxation should have higher transport levels than the experimental values to make up the shortfall, yet we did not see this.

FIG. 19.

Density and pressure profile relaxation for ions, quantities are time-averaged at the quasi steady state after saturation.

FIG. 19.

Density and pressure profile relaxation for ions, quantities are time-averaged at the quasi steady state after saturation.

Close modal

We see similar phenomenon in a different DIII-D discharge as well. Figure 20 shows the profiles of Discharge 145781 as compared to the classical shortfall case of 128913. The 145781 has higher temperature gradients at core, and we may expect the edge turbulence spreading should be less effective in increasing the core transport. This is exactly what we see in the results shown in Fig. 21. The core ITG grows quickly in the flattened edge case, and the increase of the saturation level with the edge pressure gradient drive is not big. At the outer core, Qi is still doubled by edge turbulence spreading, but Qe is almost unaffected. It is clear that the edge TEM instability has a smaller effect here.

FIG. 20.

Global simulation profiles Te, Ti, and ne of DIII-D Discharge 145781 (solid line) as compared to Discharge 128913 (dashed line).

FIG. 20.

Global simulation profiles Te, Ti, and ne of DIII-D Discharge 145781 (solid line) as compared to Discharge 128913 (dashed line).

Close modal
FIG. 21.

Simulation results for Discharge 145781. Top panel: the volume-averaged RMS δϕ; lower left panel: Qi from the global simulations compared to ONETWO; lower right panel: Qe.

FIG. 21.

Simulation results for Discharge 145781. Top panel: the volume-averaged RMS δϕ; lower left panel: Qi from the global simulations compared to ONETWO; lower right panel: Qe.

Close modal

These global results indicate that although the TEM turbulence spreading from the edge can increase the fluctuation intensity, making it nearly constant throughout the simulation domain, its effect on the transport is not as striking. Heat flux transport is still largely determined by local pressure gradients. The global effects of the edge TEM are not enough to explain the transport shortfall. Nevertheless, we note from the Discharge 145781 results that the enhancement of the core transport from edge turbulence is related to the relative strength of the edge instability. A stronger or additional edge instability that is not modeled in this study, e.g., any low-n or MHD modes, may increase the core transport more. In fact, it has been reported recently that GAMs are observed in DIII-D L-mode plasmas.38 

Only electron collisions are included in the present simulations, since the electron collision frequency is bigger than that of the ions and important for linear physics.39 Zonal flows are driven by turbulence but also regulated by nonlinear damping.40 It has been reported that ion collisions could also impact turbulence and transport via zonal flow damping.41,42 In the present case, we assume that the nonlinear damping is stronger than the collisional damping, and hence neglect the ion-ion collision. We note that collisional damping could occur over larger time scales and therefore should be considered in future shortfall studies.

To compliment recent continuum validation efforts,1,2,8,14 we have analyzed L-mode shots from both the DIII-D and Alcator C-Mod tokamaks using the nonlinear gyrokinetic particle code GEM.

For the DIII-D L-mode plasmas, the ion heat flow is comparable albeit smaller than the experimental prediction at ρ=0.5, but is much too low at location ρ=0.75. This is in agreement with earlier GYRO results.2 Linearly, both the core and edge are dominated by the ITG mode, and there is little variation of the growth rate with R/LTe, βe, and νe. By increasing the value of R/LTi by 25% and 30% at the locations ρ=0.5 and ρ=0.75, we find satisfactory agreement with the experimental heat fluxes.

For the Alcator C-Mod L-mode case, the ion heat fluxes are in satisfactory agreement with earlier GYRO results as well,8 and also within the range of experimental error. The electron heat fluxes are also in fair agreement except for the low power discharge where the electron heat fluxes are much lower than the experimental values. Linear analysis illustrates that the most unstable modes are ITG driven for the low power discharge, but both the ITG and TEM are present in the high power discharge case.

It was hypothesized in the earlier studies1 that global simulations may help reduce the discrepancy between experiment and nonlinear gyrokinetic simulation. We carry out a global simulation for the DIII-D L-mode shot with slightly modified profiles with and without the effect of edge turbulence. It appears that the core is dominated by the ITG turbulence in the absence of edge. However, when edge is included, the TEM turbulence spreading from the edge to the core enhances the ion heat flux level in the core substantially, but the overall ion heat flux is still lower than the experimentally predicted level. Therefore, we conclude that the effect of TEM edge turbulence alone is not sufficient to explain the observed shortfall in the DIII-D L-mode shots.

This material is based upon work supported in part by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences, using the DIII-D National Fusion Facility, a DOE Office of Science user facility, under Award Nos. DE-FC02-04ER54698, DOE-SC0008801, DE-FG02-07ER54917, and DE-FG02-08ER54954. DIII-D data shown in this paper can be obtained in digital format by following the links at https://fusion.gat.com/global/D3D_DMP. This research also used resources of the National Energy Research Scientific Computing Center (NERSC), which is a DOE Office of Science user facility.

We thank C. S. Chang, S. Ku, J. Lang, and R. Hager for helpful discussions. We also acknowledge fruitful discussions on “shortfall” with A. E. White and Tobias Görler.

1.
C.
Holland
,
A. E.
White
,
G. R.
McKee
,
M. W.
Shafer
,
J.
Candy
,
R. E.
Waltz
,
L.
Schmitz
, and
G. R.
Tynan
,
Phys. Plasmas
16
,
052301
(
2009
).
2.
T.
Rhodes
,
C.
Holland
,
S.
Smith
,
A.
White
,
K.
Burrell
,
J.
Candy
,
J.
DeBoo
,
E.
Doyle
,
J.
Hillesheim
,
J.
Kinsey
,
G.
McKee
,
D.
Mikkelsen
,
W.
Peebles
,
C.
Petty
,
R.
Prater
,
S.
Parker
,
Y.
Chen
,
L.
Schmitz
,
G.
Staebler
,
R.
Waltz
,
G.
Wang
,
Z.
Yan
, and
L.
Zeng
,
Nucl. Fusion
51
,
063022
(
2011
).
3.
J.
Candy
and
R. E.
Waltz
,
J. Comput. Phys.
186
,
545
(
2003
).
4.
R. E.
Waltz
,
Bull. Am. Phys. Soc.
50
,
BAPS.2012.DPP.DI3.2
(
2012
).
5.
R. E.
Waltz
,
Phys. Plasmas
20
,
012507
(
2013
).
6.
W. M.
Stacey
,
Phys. Plasmas
20
,
112503
(
2013
).
7.
A. E.
White
,
N. T.
Howard
,
M.
Greenwald
,
M. L.
Reinke
,
C.
Sung
,
S.
Baek
,
M.
Barnes
,
J.
Candy
,
A.
Dominguez
,
D.
Ernst
,
C.
Gao
,
A. E.
Hubbard
,
J. W.
Hughes
,
Y.
Lin
,
D.
Mikkelsen
,
F.
Parra
,
M.
Porkolab
,
J. E.
Rice
,
J.
Walk
,
S. J.
Wukitch
, and
Alcator C-Mod Team
,
Phys. Plasmas
20
,
056106
(
2013
).
8.
N.
Howard
,
A.
White
,
M.
Reinke
,
M.
Greenwald
,
C.
Holland
,
J.
Candy
, and
J.
Walk
,
Nucl. Fusion
53
,
123011
(
2013
).
9.
N. T.
Howard
,
A. E.
White
,
M.
Greenwald
,
M. L.
Reinke
,
J.
Walk
,
C.
Holland
,
J.
Candy
, and
T.
Görler
,
Phys. Plasmas
20
,
032510
(
2013
).
10.
N. T.
Howard
,
A. E.
White
,
M.
Greenwald
,
C.
Holland
, and
J.
Candy
,
Phys. Plasmas
21
,
032308
(
2014
).
11.
N.
Howard
,
C.
Holland
,
A.
White
,
M.
Greenwald
, and
J.
Candy
, “Synergistic cross-scale coupling of turbulence in a tokamak plasma,”
Phys. Plasmas
(submitted,
2014
).
12.
F.
Jenko
,
W.
Dorland
,
M.
Kotschenreuther
, and
B. N.
Rogers
,
Phys. Plasmas
7
,
1904
(
2000
).
13.
T.
Görler
,
X.
Lapillonne
,
S.
Brunner
,
T.
Dannert
,
F.
Jenko
,
F.
Merz
, and
D.
Told
,
J. Comput. Phys.
230
,
7053
(
2011
).
14.
D.
Told
,
F.
Jenko
,
T.
Görler
,
F. J.
Casson
,
E.
Fable
, and
ASDEX Upgrade Team
,
Phys. Plasmas
20
,
122312
(
2013
).
15.
T.
Görler
, private communications (2014).
16.
Y.
Chen
and
S. E.
Parker
,
J. Comput. Phys.
189
,
463
(
2003
).
17.
Y.
Chen
and
S. E.
Parker
,
J. Comput. Phys.
220
,
839
(
2007
).
18.
R. V.
Bravenec
,
Y.
Chen
,
J.
Candy
,
W.
Wan
, and
S.
Parker
,
Phys. Plasmas
20
,
104506
(
2013
).
19.
Y.
Chen
,
S. E.
Parker
,
W.
Wan
, and
R.
Bravenec
,
Phys. Plasmas
20
,
092511
(
2013
).
20.
Z.
Lin
and
T. S.
Hahm
,
Phys. Plasmas
11
,
1099
(
2004
).
21.
T. S.
Hahm
,
P. H.
Diamond
,
Z.
Lin
,
K.
Itoh
, and
S.-I.
Itoh
,
Plasma Phys. Controlled Fusion
46
,
A323
(
2004
).
22.
T. S.
Hahm
,
P. H.
Diamond
,
Z.
Lin
,
G.
Rewoldt
,
O.
Gurcan
, and
S.
Ethier
,
Phys. Plasmas
12
,
090903
(
2005
).
23.
W. X.
Wang
,
Z.
Lin
,
W. M.
Tang
,
W. W.
Lee
,
S.
Ethier
,
J. L. V.
Lewandowski
,
G.
Rewoldt
,
T. S.
Hahm
, and
J.
Manickam
,
Phys. Plasmas
13
,
092505
(
2006
).
24.
S.
Ku
,
J.
Abiteboul
,
P.
Diamond
,
G.
Dif-Pradalier
,
J.
Kwon
,
Y.
Sarazin
,
T.
Hahm
,
X.
Garbet
,
C.
Chang
,
G.
Latu
,
E. S.
Yoon
,
P.
Ghendrih
,
S.
Yi
,
A.
Strugarek
,
W.
Solomon
, and
V.
Grandgirard
,
Nucl. Fusion
52
,
063013
(
2012
).
25.
P.
Diamond
,
S.
Champeaux
,
M.
Malkov
,
A.
Das
,
I.
Gruzinov
,
M.
Rosenbluth
,
C.
Holland
,
B.
Wecht
,
A.
Smolyakov
,
F.
Hinton
,
Z.
Lin
, and
T.
Hahm
,
Nucl. Fusion
41
,
1067
(
2001
).
26.
W.
Wan
,
S. E.
Parker
,
Y.
Chen
, and
F. W.
Perkins
,
Phys. Plasmas
17
,
040701
(
2010
).
27.
W.
Wan
,
S. E.
Parker
,
Y.
Chen
,
G.-Y.
Park
,
C.-S.
Chang
, and
D.
Stotler
,
Phys. Plasmas
18
,
056116
(
2011
).
28.
D.
Smith
,
S.
Parker
,
W.
Wan
,
Y.
Chen
,
A.
Diallo
,
B.
Dudson
,
R.
Fonck
,
W.
Guttenfelder
,
G.
McKee
,
S.
Kaye
,
D.
Thompson
,
R.
Bell
,
B.
LeBlanc
, and
M.
Podesta
,
Nucl. Fusion
53
,
113029
(
2013
).
29.
Y.
Chen
,
S. E.
Parker
,
G.
Rewoldt
,
S.-H.
Ku
,
G.-Y.
Park
, and
C.-S.
Chang
,
Phys. Plasmas
15
,
055905
(
2008
).
30.
W.
Wan
,
S. E.
Parker
,
Y.
Chen
,
Z.
Yan
,
R. J.
Groebner
, and
P. B.
Snyder
,
Phys. Rev. Lett.
109
,
185004
(
2012
).
31.
W.
Wan
,
S. E.
Parker
,
Y.
Chen
,
R. J.
Groebner
,
Z.
Yan
,
A. Y.
Pankin
, and
S. E.
Kruger
,
Phys. Plasmas
20
,
055902
(
2013
).
32.
A. E.
White
,
L.
Schmitz
,
G. R.
McKee
,
C.
Holland
,
W. A.
Peebles
,
T. A.
Carter
,
M. W.
Shafer
,
M. E.
Austin
,
K. H.
Burrell
,
J.
Candy
,
J. C.
DeBoo
,
E. J.
Doyle
,
M. A.
Makowski
,
R.
Prater
,
T. L.
Rhodes
,
G. M.
Staebler
,
G. R.
Tynan
,
R. E.
Waltz
, and
G.
Wang
,
Phys. Plasmas
15
,
056116
(
2008
).
33.
I.
Manuilskiy
and
W.
Lee
,
Phys. Plasmas
7
,
1381
(
2000
).
34.
Y.
Chen
and
S.
Parker
,
Phys. Plasmas
8
,
2095
(
2001
).
35.
R. V.
Bravenec
,
J.
Candy
,
M.
Barnes
, and
C.
Holland
,
Phys. Plasmas
18
,
122505
(
2011
).
36.
See http://w3.pppl.gov/transp/ for more information about TRANSP.
37.
R. L.
Miller
,
M. S.
Chu
,
J. M.
Greene
,
Y. R.
Lin-Liu
, and
R. E.
Waltz
,
Phys. Plasmas
5
,
973
(
1998
).
38.
G.
Wang
,
W. A.
Peebles
,
T. L.
Rhodes
,
M. E.
Austin
,
Z.
Yan
,
G. R.
McKee
,
R. J. L.
Haye
,
K. H.
Burrell
,
E. J.
Doyle
,
J. C.
Hillesheim
,
M. J.
Lanctot
,
R.
Nazikian
,
C. C.
Petty
,
L.
Schmitz
,
S.
Smith
,
E. J.
Strait
,
M. V.
Zeeland
, and
L.
Zeng
,
Phys. Plasmas
20
,
092501
(
2013
).
39.
M.
Kotschenreuther
,
G.
Rewoldt
, and
W.
Tang
,
Comput. Phys. Commun.
88
,
128
(
1995
).
40.
M. N.
Rosenbluth
and
F. L.
Hinton
,
Phys. Rev. Lett.
80
,
724
(
1998
).
41.
Z.
Lin
,
T. S.
Hahm
,
W. W.
Lee
,
W. M.
Tang
, and
P. H.
Diamond
,
Phys. Rev. Lett.
83
,
3645
(
1999
).
42.
Y.
Xiao
,
P. J.
Catto
, and
K.
Molvig
,
Phys. Plasmas
14
,
032302
(
2007
).