Proceedings of ITC/ISHW2007

Development of non-local neoclassical transport code for helicalconfigurations
Shinsuke Satake1) , Ryutaro Kanno1,2) and Hideo Sugama1,2)
1) National Institute for Fusion Science,322-6 Oroshi-cho, Toki 509-5292,, Japan
2) Department of Fusion Science, Graduate University for Advanced Studies, Toki 509-5292, Japan

The progress in 3-dimensional, non-local neoclassical transport simulation code “FORTEC-3D” is described.
The main purpose of the code is to solve the drift-kinetic equation in general 3-dimensional configuration using
the δ f Monte-Carlo method, and to calculate neoclassical fluxes and the time evolution of the ambipolar radial
electric field simultaneously. In this article, we explain the new numerical schemes adopted for FORTEC-3D inorder to overcome problems happened especially in the cases where the bifurcation of radial electric field occurs.
Examples of test simulation for a LHD magnetic field configuration with bifurcated electric field are also shown.
With improved numerical schemes, FORTEC-3D can calculate neoclassical fluxes and trace the time evolution
stably as long as several ion collision times, which is long enough to observe GAM damping and transition of
the ambipolar electric field.
Keywords: neoclassical transport, δ f Monte-Carlo method, ambipolar radial electric field

1Introduction

Detailed calculation of neoclassical transport in 3dimensional configuration plasmas such as LHD is important for transport analysis, since in a 3-dimensional system the ambipolar condition Γi (ρ , E) = Γe (ρ , E) determines the radial electric field profile E(ρ ), where Γi andΓe are ion and electron particle fluxes across the flux surface ρ = const. It is known that the ambipolar condition
is mainly determined by neoclassical transport. If the ambipolar condition has multiple roots, bifurcation of radial
electric field profile occurs.[1] In LHD plasmas, appearance of positive electric field, or the “electron root”, ispreferable since it reduces neoclassical transport comparedwith that in negative one, or the “ion root”.[2] Moreover,strong sheared electric field profile at the bifurcation pointis generally considered to be favorable from the viewpointof suppression of anomalous transport.
Neoclassical transport theory for 3-dimensional helical configuration[3, 4] has been constructed under the assumption of local transport model where the typical orbit width in the minor-radius direction is assumed to benegligible compared with the background gradient scale
lengths of a plasma. Recently non-local effects, or finiteorbit-width (FOW) effects, on neoclassical transport haveattracted attention in the analysis of core transport in tokamaks, where particle orbits with large widths break the assumption of conventional neoclassical theory.[5, 6] In helical configurations, the deeply ripple-trapped and transition particles show large deviation in radial direction. TheFOW effect of these orbits will be important in neoclassical transport analysis if the collisionality becomes quite
author’s e-mail: satake@nifs.ac.jp

low. The presence of strong-sheared electric field will also
brings the non-local effect to neoclassical transport. However, the applicability of analytic neoclassical theory tothese situations is questionable since it neglects the FOW
effect from the beginning.
From the considerations above, we have been developing a simulation code to solve the drift-kinetic equation(DKE) including the FOW effect in 3-dimensional configuration. The simulation code, FORTEC-3D,[7, 8] uses theδ f Monte-Carlo method[9, 10] which has been applied tosome other transport codes both for tokamaks[11] and forhelical configurations.[12, 13] The features of FORTEC3D are : (1) It uses a conserved-form linearized FokkerPlanck collision operator. (2) It is a global simulation. The
entire confinement region is solved at once. (3) Time evolution of radial electric field is solved simultaneously. The
ambipolar electric field is then determined in consistency
with neoclassical fluxes. (4) To reduce calculation time,
only the ion transport is solved by the δ f method, andelectron transport is solved by GSRAKE code[14] whichsolves bounce-averaged DKE. Note that GSRAKE solution does not include the FOW effect, and then only thenon-local effect for ions transport is treated in FORTEC3D.
So far, FORTEC-3D has successfully applied for LHD
configurations to solve the formation of ambipolar field
for ion roots,[7] and to study the configuration dependenceof GAM oscillation and damping[8]. However, from several test calculations we found that it was difficult to apply
FORTEC-3D to the cases with electron roots, since bifurcated radial electric field profile became unphysical shape
as shown later. We also found that numerical noise in particle flux and numerical error in collision operator were in-Proceedings of ITC/ISHW2007
tolerably large, which had not found in tokamak cases. We
ascertained the cause of these problems were from inaccuracy in radial grids for evaluating flux and electric field,
and from the emergence of huge-weight markers.
In this paper, improvements for numerical schemes inFORTEC-3D recently applied to overcome the problems
above are explained. In section 2, basic equations for δ fMonte-Carlo method are reviewed. In section 3, modification to collision operator is described. The new operator has good conservation property with less marker numbers. In section 4, improvement in the evaluation of fluxand electric field are explained. Adoption of staggeredmesh in radial direction to evaluate these two quantities
makes it possible to simulate the formation of bifurcatedradial electric field profile stably. In section 5, new filtering
scheme for marker weights to reduce numerical noise is introduced. By comparing several tests with varied strengthof filters, it is shown that the filtering scheme does not affect the solutions. By these improvements, now FORTEC3D is ready to solve neoclassical transport in helical plasmas in varied profiles and simulate evolution of electricfield including bifurcations.

2 Basic equations of the δ f method
In the δ f method, time development of the perturbationof plasma distribution function from the local Maxwellianδ f = f − fM is solved according to the DKE´Dδ f∂δ f ³≡+ vk + vd · ∇δ f −Ct p (δ f )Dt∂t= −vd · ∇ fM + P fM ,(1)where Ct p and P are test-particle and field-particle partsof linearized collision operator. The magnetic field isgiven in the Boozer coordinate system (ψ , θ , ζ )[15] as
B = ∇ψ × ∇θ + ι ∇ζ × ∇ψ . Practically, we use ρ for anormalized
p radial coordinate defined from toroidal flux ψas ρ = ψ /ψout . A MHD equilibrium magnetic field isconstructed from the VMEC code[16]. The time evolutionof radial electric field E = −dΦ/d ρ ∇ρ = Eρ ∇ρ is solved
from the following equation

ε0 ε⊥

∂ Eρ
= −e [zi Γi − Γe ] ,∂t(2)wheresubscripts
i and® e describe
particle species, and ε⊥ ≡£­® ­¤|∇ρ |2 + c2 |∇ρ |2 /v2A .
To solve eq. (1), two weights, w and p, are introduced
which satisfy the relations wg = δ f and pg = fM , respectively, where g is 5-dimensional simulation marker distribution function. Since the time evolution of marker distribution is given by Dg/Dt = 0, we obtaindw
dt
dpdt

==

¤p £−vd · ∇ + P fM ,fM
pv · ∇ fM .
fM d(3a)(3b)

The numerical procedures for the collision operator and foreq. (2) are described in the following sections.

3Collision operator

The linearized Fokker-Planck operator in FORTEC-3D is
made to satisfy the following relations,Z¡¢d 3 vM Ct p (δ f ) + P fM = 0 for M = {1, vk , v2 }, (4)Ct p (δ f ) + P fM = 0 for δ f = (c0 + c1 · v + c2 v2 ) fM . (5)The test-particle operator Ct p is expressed by random scattering in the (vk , v⊥ )-space. The field-particle operator isgiven as followsr¸·π
P = a 1−3(φ − φ 0 )x−1/2 δ n
2+ bvk x−3/2 φ δ P + cx−1/2 (φ − φ 0 )δ E,2 , φ (x) is the error function, andwhere x ≡ v2 /vthZno{δ n, δ P, δ E} = d 3 v 1, vk , v2 Ct p (δ f )

(6)

(7)

are changes in constants-of-motions by Ct p only. Previously, the relation (3b) has been used to determine constants (a, b, c) = (1/n, 2/nvth , 2/3nvth ). However, it was
found that the numerical error in the conservation low (4)became larger as the numerical noise fM − pg becamelarger. The numerical error is significant in 3-D helicalconfiguration cases compared with 2-D cases, because in3-D cases it is difficult to keep enough marker population
in a unit volume as in 2-D cases, and the marker distribution g is distorted in the velocity space in the presenceof ripple-trapped particles. In the improved version, wedo not use (3b) any more to make P fM , but the factors(a, b, c) are determined at every moment of collision sothat the conservation low is strictly satisfied. In fig. 1,
the residual relative error in momentum and energy arisen
at one operation of Ct p + P fM are compared between oldand new P fM for several tests with varied marker numbers. For previous operator, the numerical error is larger
for fewer marker case. However, the error in new operator
is just the rounding-error level even in the 1600-markers
calculation. We have also confirmed the other property ofcollision operator (5) is also kept in the new operator.

4 Radial meshes
The relation between neoclassical fluxes and radial electricfield (2) is solved on discrete meshes in the ρ -coordinate.
Γi is evaluated by the volume averaged value between i-thand i + 1-th meshes as follows
ZZ
1d 3 x d 3 vρ̇δ f
Γi (ρ 1 ) =i+ 2∆Vi+1/2 ∆Vi+1/2=1∆Vi+1/2 {k|ρ ≤∑ρ <ρiki+1 }wk ρ̇kC(ρk , i),(8)Proceedings of ITC/ISHW2007  !
"$#&% ')(&*,+ - -6

Residual err in energy (log)-8n=64000,oldn=16000,oldn=4000,oldn=1600,oldn=64000,new
n=1600,new







-10



-12




-14

    

ρ

-16
-16

-14

-12-10-8
Residual err in moment (log)-6-4Fig. 1 Relative error in momentum and energy in one operation of Ct p + P fM . The time interval for the operationis taken ∆t = 2 × 10−3 τii .

where k is the marker index, and C(ρk , i) is a shaping factor
of a marker.[17] For electron flux GSRAKE solution is referred to : Γe (ρi+1/2 ) = Γe (ρi+1/2 , Eρ (ρi+1/2 ,t))GSRAKE .
Previously, the electric field and fluxes were evaluated on
the same points ρ = ρi+1/2 according to eq. (2). It is illustrated in Fig. 2. However, using this scheme, we foundthat the electric field profile fell into an unphysical shape asshown in Fig. 3, where multiple roots for ambipolar condition are expected from GSRAKE solutions. This is because all three values in eq. (2) were evaluated on the sameposition, that is, the evaluation was relied on a local information too much. Therefore, we adopted a new scheme
which uses staggered-mesh for Eρ and Γi , Γe as illustrated
in Fig. 2. Fluxes are evaluated on half-grids ρ = ρi±1/2
and time evolution of Eρ (ρi ) is calculated as follows

ε0 ε⊥¤£∂ Eρ (ρi )= −e zi Γ̄i (ρi ) − Γ¯e (ρi ) ,∂t(9)where Γ̄i (ρi ) and Γ¯e (ρi ) are averaged values of those which
evaluated at ρ = ρi±1/2 . With this new scheme, FORTEC3D can simulate continuous transition from ion-root toelectron-root as shown in Fig. 3.
(new)
(old)

Fig. 2 Radial meshes and positions on which fluxes and electricfield are evaluated in old and new schemes.

Fig. 3 Comparison of electric field profile between old and newschemes for flux and electric field. Squares are estimationof ambipolar field from GSRAKE.

5 Filtration
Reducing the numerical noise without relying on a massive
number of markers is essential to any Monte-Carlo simulation. In the δ f method, the origin of big noise is frommarkers which has huge weights w or p. If collisions andelectric potential are neglected, the marker weights can bedetermined as
µp1
w1

¶3/2

T0ne
p0 1
n0 T1
= w0 − (p1 − p0 ),=" µ¶2 µ 2 ¶2 #mv
mv2− 2T1 + 2T0
10, (10a)(10b)where subscript 0 and 1 represents initial and present values of them at marker position. From eq.(10), it is expectedthat markers which travels long distance in radial directionwill have large marker weights. It is inevitable to appear
such large-drift particles in low-collisionality helical plasmas, we need to set limits for weight values. In reality,
we found only less than 0.1% of simulation markers havehuge weights p1 ∼ |w1 | À 10p0 , and these markers makes
the calculation very noisy and unstable. Therefore, we apply filters for markers by setting limits for p, |w|/p, andv/vth . The second limiter is following from the assumption in δ f method that δ f / fM ¿ 1, and the third limiteris because fast ions tend to have large orbit widths. Themarkers which breaks the limits are filtrated out and redistributed around the magnetic axis, which contributes tokeep marker population there.
To check the effect of filters, we have carried out several tests with varied strength of filters as shown in Table
1. Note here that the “old” simulation used old P operator and was not equipped with filters for p. In fig. 4,
comparison of radial electric field at three different times
with varied filters are shown. As is shown there, the timeevolution of Eρ and quasi-steady state profiles (at the lastfigure of three) are almost the same regardless of strengthof filters. To check the time evolution in detail, we compare Eρ and Γi on ρ = 0.45 for three filters in Fig. 5. Here,error-bars are evaluated from the variance between every
30 steps. It is found that the filters successfully suppressProceedings of ITC/ISHW2007
the numerical noise even in the weakest filter case, though
the start timing of transition differs slightly among threesimulations. It is also notable that the simulation marker
number can be reduced if we adopt the filters.

Fig. 5 Comparison of time evolution of Er and Γi between old,filter 1 and filter 5 calculations on ρ = 0.45 surface. The
oscillation at the beginning phase is GAM.

Table 1 Filter strengths and marker numbers used in test-runs.

oldfilter 1
filter 2
filter 3
filter 4
filter 5Fig. 4 Comparison of electric field profile between old and newschemes for different filters.

6 Conclusion
The improvements for numerical schemes in FORTEC-3D
were proved to reduce numerical errors and noises significantly with little changes in the observable values such as
flux and electric field, etc. FORTEC-3D will be applied to
study the FOW effects in helical plasmas, bifurcation phenomena, and so on. We are also planning to apply the δ fmethod to solve electron transport and evaluate the bootstrap current in the presence of ambipolar electric field.
[1] D. E. Hastings et al.,Nuclear Fusion 25, 445 (1985).
[2] R. Kanno et al.,et al.,Nuclear Fusion 37, 1463 (1997).

p58102510|w|/p
2.0
2.0
2.5
2.54.010.0v/vth
5.0
4.5
4.5
4.5
4.5
4.5

marker num.
3.84 × 107
3.84 × 107
3.84 × 107
1.28 × 107
1.92 × 107
1.92 × 107[3] N. Nakajima et al.,J. Phys. Soc. Japan 61, 833 (1992).
[4] M. Wakatani, Stellarators and Heliotron Devices (Oxford
Univ. Press, 1998), Chap. 7.
[5] P. Helander, Physics of Plasmas 7, 2878 (2000).
[6] S. Satake et al.,Physics of Plasmas 9, 3946 (2002) .
[7] S. Satake et al.,Plasma Fusion Res. 1, 002 (2006).
[8] S. Satake et al.,Nuclear Fusion 47, 1258 (2007).
[9] W. X. Wang et al.,Plasma Phys. Control. Fusion 41 1091
(1999).
[10] S. Brunner et al.,Physics of Plasmas 6, 4504 (1999).
[11] W. X. Wang et al.,Physics of Plasmas 13 082501 (2006).
[12] J. L. Lewandowski et al.,Physics of Plasmas 8, 2849 (2001).
[13] M. Yu. Isaev et al.,Fusion Science and Technology 50, 440(2006).
[14] C. D. Beidler et al.,Plasma Phys. Control. Fusion 43, 1131
(2001).
[15] A. H. Boozer, Phys. Fluids 23, 904 (1980).
[16] S. P. Hirshman et al.,J. Comput. Phys. 96, 99 (1991).
[17] C. K. Birdsall et al.,Plasma Physics via Computer Simulation (Inst. Phys. Pub., 1991), Chap. 4.