J. Plasma
J.
Plasma Fusion
Fusion Res.
Res. SERIES,SERIES,Vol.
Vol.6 6(2004)(2004)203–206
000–000

Development of a Neoclassical Transport Database
by Neural Network Fitting in LHD
WAKASA Arimitsu, MURAKAMI Sadayoshi1, YAMADA Hiroshi2, YOKOYAMA Masayuki2,MAASSBERG Hening3, BEIDLER Craig D.3, WATANABE Kiyomasa2,
OIKAWA Shun-ichi and LHD Experimental Group2
Graduate School of Engineering, Hokkaido University, Sapporo, 060-8628, Japan1
Department of Nuclear Engineering, Kyoto University, Kyoto, 606-8501, Japan
2
National Institute for Fusion Science, Toki, 509-5292, Japan3Max-Planck-Institute für Plasmaphysik, EURATOM Ass., Greifswald, Germany(Received: 18 December 2003 / Accepted: 19 May 2004)

Abstract
A database of neoclassical transport coefficients for the Large Helical Device is developed using normalized
mono-energetic diffusion coefficients evaluated by Monte Carlo simulation code; DCOM. A neural network fitting
method is applied to take energy convolutions with the given distribution function, e.g. Maxwellian. The databasegives the diffusion coefficients as a function of the collision frequency, the radial electric field and the minor radius
position.
Keywords:
LHD, neoclasiscal transport, neural network, Monte Carlo method

1. Introduction

we have developed the DCOM (Diffusion COeficient
calculator by Monte carlo method) code [8].
The calculation of the mono-energetic diffusion
coefficient using the Monte Carlo method, however, takes
very long time. Since the neoclassical diffusion coefficientdepends on the radial position, the radial electric field, andthe collision frequency, it is difficult to calculate wheneverthose parameters change.
All the elements of the transport matrix can be obtained
from the convolution of Maxwellian distribution function withthree different mono-energetic coefficients [9];Because of the radial drift motions of helical trapped
particles, the neoclassical transport increases as the collision
frequency decreases (1/ν regime) in the long-mean-free-path(LMFP) regime in heliotrons. The neoclassical transport
would play an important role as well as the anomaloustransport in fusion plasma, where the confinement of hightemperature plasma is required. Also, recent LHD experimentresults suggest a key role of the neoclassical transport in theelectron internal transport barrier [1] formation to determineradial electric field.
The estimate of the neoclassical transport coefficient is
an important issue, and many studies have been done toevaluate the neoclassical transport coefficient analytically andnumerically in helical systems. Among them the DKES (Drift
Kinetic Equation Solver) code [2,3] has been commonly used
for the experimental data analysis [4,5] and for theoreticalpredictions [6,7]. However, in the LMFP regime, a large
number of Fourier modes must be used for the distributionfunction and a convergence problem can be seen.
Monte Carlo method is also used to evaluate theneoclassical transport coefficient, where the diffusion
coefficients are estimated by the radial diffusion of test
particles. This method does not have a convergence problem
and is more applicable in the LMFP regime. Owing to above,kj

D =4π∫2j⎡ ⎛ v ⎞ 2 ⎤ dv
⎛ v ⎞⎢⎥,D (v) ⎜⎟ exp ⎢ − ⎜⎟ ⎥vvvth⎝ th ⎠⎝
th ⎠⎣
⎦k3

⎛ v ⎞D (v) = D D* (ν *, G ) ⎜⎟ ,⎝ vth ⎠kkP

where k = e, i (for electrons or ions) and j = 1, 2, 3 (relatedto the diffusion, the bootstrap current and parallel resistivity),
DP is the tokamak plateau value of mono-energetic case, DP= (π /16)(υ3/ιRωc2), R is the major magnetic radius, υ is thevelocity of test particles, ι is the rotational transform
normalized by 2π, ωc is cyclotron frequency of test particles,

Corresponding author’s e-mail: wakasa@fusion.qe.eng.hokudai.ac.jp1203©2004 by The Japan Society of Plasma
Science and Nuclear Fusion ResearchWakasa A. et al., Development of a Neoclassical Transport Database by Neural Network Fitting in LHD

ν* is the collision frequency normalized by vι/R, vth is thethermal velocity, D* is the diffusion mono-energetic
coefficient normalized by Dp , G is the normalized radial
electric field, G = (R/ιr)(Er /υB), which measures the effectof the poloidal rotation induced by the radial electric field tothe poloidal rotation of circulating particles, r is the minor
radial position, Er is the radial electric field and B is themagnetic field strength. In order to calculate the transportcoefficient for Maxwellian distribution, it is necessary tointerpolate the mono-energetic diffusion coefficient calculatedby DCOM code. A conventional fitting method assuming thecommon analytical relations depending on each collisionality
regimes with their matching conditions, has been used so far.
However, this technique can’t describe accurately theoverlapping of the each collisionality regimes.
Recently, the technique of neural networks (NNW) [10],that has strong nonlinear and has high fitting abilities havebecome a center of attraction. The NNW has their origins instudies of the brain, and has been the subject of active
research since the mid 1940’s. In the last few years inparticular, the NNW has been used to solve a large range ofpractical problems in many areas of physics and engineering.
It is appropriate to use the NNW for the interpolation of thediffusion coefficient in LHD.
In this paper, we show the neoclassical transport databasefor LHD using NNW fitting method. The evaluation of themono-energetic diffusion coefficients by DCOM code isexplained for the LHD standard configuration (Rax = 3.75 m)
in Sec. 2. In Sec. 3, we have constructed the data basecontaining the normalized mono-energetic diffusion
coefficients for several radial points, radial electric field andcollisionalities by using the NNW.

4102

D*

100

10-2

10-4

10Fig. 1-710-5

10-3

10-1

110
ν*103

10The monoenergetic diffusion coefficient, D* vs thenormalized collision frequency, ν * in standard LHD
configuration, R ax = 3.75 m for different values ofelectric potential corresponding to G = 0, 1 × 10-3, 3 ×10-3, 1 × 10-2, 3 × 10-2, 1 × 10-1 and 3 × 10-1 (crosses,squares, circles, triangles, down triangles, diamondsand pentagons, respectively; closed signs arecalculating by DCOM and opened signs are by DKES).

times, τ, the diffusion coefficient can be evaluated by taking
the mean square displacement of N particles as

D=1∑ ⎛⎝ri − r ⎞⎠ ,2Nτ
N

2i =1

where ri stands for the radial position of i-th particle, andr = ( 1/ N ) ∑ i =1 ri .
N

2. Calculation of the monoenergetic
coefficients by Monte Carlo method

In order to obtain mono-energetic diffusion coefficients,
we have chosen N = 1000 with the energy of 1.0 × 10-3 eV.
The magnetic field is set to be 3 T at the magnetic axis. The
test particles are followed for several collisional times until
the evaluated diffusion coefficient is converged.
In this work, we calculated the normalized monoenergetic diffusion coefficients; standard LHD configurations,
Rax = 3.75 m, 5 radial positions normalized the minor radius,ρ = 0.1, 0.25, 0.5, 0.75, 0.9, 14 normalized collisionalities,
3.16 × 10-2 < ν* < 1.00 × 106 and 7 electric fields, 0.00 < G< 0.03.
The calculation results obtained by DCOM code without
E r are shown by the diagonal cross sign in Fig. 1. The
normalized diffusion coefficients D* as a function of ν* at ρ= 0.5 in the standard configuration are shown. Shown in Fig.
1 are also the DKES results (cross sign) and analytical results(dashed line). The results of DCOM and DKES agree wellfrom Pfirsch-Schlüter (P-S) regime to 1/ν regime for bothconfiguration cases. The dashed lines represent analyticalresult of diffusion coefficient Dana (= Da + D1/ν ) which is thesum of non-axisymmetric contribution D1/ν given by a multihelicity model [11] and the axisymmetric contribution Da

The neoclassical transport coefficient has been evaluatedusing Monte Carlo technique directly following the particle
orbits. We evaluate a mono-energetic local diffusion
coefficient, D. In the simulation the N monoenergetic particles
are released at an initial minor radial position, r0 , where theparticles are randomly distributed in the poloidal, toroidal
coordinates, and in the pitch angle space. The test particle
orbits are followed solving the guiding-center equations in theBoozer coordinates with 50 Fourier modes of the magneticfield. The magnetic field in the Boozer coordinates isconstructed from the magnetohydrodynamics (MHD)
equilibrium obtained by the VMEC code. The pitch angle
scatterings are taken into account in the Monte Carlo collision
operator based on the binominal distribution [10]. The pitch
angle scattering after the time interval ∆t is given in terms ofλ(= υ || /υ) by1
22
λ n + 1 = λ n (1 − ν ∆ t ) + σ ⎡ ⎛⎝1 − λ n ⎞⎠ ν ∆ t ⎤ ,⎣⎦where σ takes +1 or –1 with equal probabilities, ν is thedeflection collision frequency. After several characteristic2204

Wakasa A. et al., Development of a Neoclassical Transport Database by Neural Network Fitting in LHD

given byoutput
3/2
3/2
Da = ⎛⎝ Dbp + DPS ⎞⎠D*

2/3,

ω2

whereDbp = Db Dp⎛D + D ⎞ ,p⎠⎝ bhidden
units

1/2
22
Db = ⎛⎝B1, 0 ε t ⎞⎠ ⎛⎝vd R ωc ι ⎞⎠ ν ,bias

ω

Dp = ⎛⎝ 2 5⎞⎠ ⎛⎝B1, 0 ε t⎞⎠ ⎛⎝vd v ωcι⎞⎠1bias
andν*inputs

DPS = ⎛⎝ 7 5⎞⎠ ⎛⎝ B1, 0 ε t⎞⎠ ⎛⎝ vd R ωc ι ⎞⎠ ν ,Gρ

2Fig. 2where υd = υ2/RB2 is the drift velocity of the test particles.
As shown Fig. 1, the correct diffusion coefficient using
analytical results in plateau regime can not be obtained.
Next, we considered the effect of the radial electric field,Er . In Fig. 1, the diffusion coefficients at ρ = 0.5 are shownfor different amplitude of radial electric field. The open signs
represent the computational results by DKES and the closed
signs represent by DCOM. The strong reductions in thediffusion coefficients as the radial electric field strengthens
are observed in the 1/ν regime.

Schematic view of the MLP1 neural network used in thecalculations of the normalized D*.

Averaged relative error [%]100

3. Construction of the database of thediffusion coefficient using Neural Net
Work

10

10The diffusion coefficient of LHD has a strong
nonlinearity and is a function of several parameters; ν*, Er ,ρ, Rax , β, and etc. The NNW is suitable for the interpolationof the nonlinear function which depends on the severalparameters. We consider the most widely used NNW, known
as a multilayer perceptron (MLP), with only one hidden layer,
MLP1 to construct the database of the neoclassical transport.
The outputs of this network yn(n = 1, 2, ...N ) can be written
as a function of its inputs xm(m = 1, 2, ...M) through a set of
coefficients, or weights, {ω 1hm , b1h , ω 2hn , b2n}(h = 1,2, ...H), asFig. 3

5

10 15 20 25 30 35 40
number of hidden unit of NNW45

A plot of the averaged relative error versus the numberof the hidden unit.

the underlying trend in the data and which therefore gives
poor predictions for a new data. Figure 3 shows a plot of theaveraged relative error versus the order of the hidden units.
The results in Fig. 3 decided a number of the hidden unit
was 14.
We can obtain D* by using the NNW for arbitrary ν*, Gand ρ. The NNW results at ρ = 0.5 for Rax = 3.75 m areshown in Fig. 4. The square symbols represent thecomputational results by DCOM with the training of theNNW and the solid lines represent the outputs from NNW.
The NNW results agree well with DCOM results in eachregion of 1/ν, plateau, and P-S, and smooth fitting is done in
the joint of each region. The triangle symbols show theDCOM results without the training of NNW and the dashedlines show the results with NNW. It should be noted that theaccuracy of the NNW interpolation with respect to the electricfield is also good.
Figure 5(a) shows contour plot of D*, obtained by theNNW, at ρ = 0.5. The horizontal axis indicates ν* and thevertical axis is G. We found that D* does not depend on theelectric field in P-S regime and D* is improved in 1/ν region⎛ H 2 ⎛ M 11⎞2⎞yn (x1, x 2 , ..., x M ) = f ⎜ ∑ ωhn f ⎜ ∑ ωhm x m + b n ⎟ + bn ⎟⎠⎝ h =1⎠⎝ m =1
with a sigmoid function f(x) ≡ tanh(x). In this work, the NNW
was constructed as in Fig. 2. The inputs are ν*, G and ρ andthe output is D*. To equate the computational result of
DCOM and output of the NNW, the weights ω have beencorrected. The BFGS (Broyden-Fletcher-Goldfarb-Shanno)
method [12] of the quasi-Newton method was used for thisminimum value search in this study.
The accuracy of the NNW depends on the number ofthe hidden unit. The error margin decreases as a number ofthe hidden unit increases. However, if the order of the hidden
unit is increased too far, the phenomenon of the overfitting
occurs which gives a very small error with respect to thetraining data, but which again gives a poor representation of3205

Wakasa A. et al., Development of a Neoclassical Transport Database by Neural Network Fitting in LHD

as the electric field strengthens. Figure. 5(b) shows thecontour plot of D* as a function of ν* (horizontal axis) and ρ(vertical axis). It is seen that the plateau regime has narrowed
as the radial position goes toward outside of the plasma ineach configuration.

necessary to interpolate the values of D* as a function of thecollision frequency, the radial electric field, the minor radiusand etc. The neural network (NNW) is applied to fit thediffusion coefficient of LHD, which shows complex behavior
in the several collisional regimes; i.e. ν, ⎯ν√ , 1/ν, plateau, andP-S regimes. A multilayer perceptron NNW with only one
hidden layer, usually known as MLP1, was used. Input
parameters are ρ, ν* and G, and D* can be obtained as anoutput of the NNW.
A single set of D* database can be applicable to theevaluation of the transport coefficients of one configurationin any density and temperature parameter plasma. Also, theNNW database can evaluate the mono-energetic diffusion
coefficient in a very short time and we can evaluate theneoclassical transport including the ambipolar radial electricfield in a minute using a present PC. On the other hand, theconfiguration changes significantly in a finite beta plasma ofLHD and, thus, the construction of D* database in the finite
beta is necessary for accurate evaluation of neoclassicaltransport. The database construction in a finite beta is inprogress.

4. SummaryD*

We have constructed the database of the neoclassicaltransport in LHD. The mono-energetic neoclassical diffusion
coefficients are evaluated numerically using Monte Carlotechnique. We have used a Monte Carlo simulation code,DCOM, in which the neoclassical diffusion coefficients areestimated by the radial diffusion of test particles in Boozer
co-ordinates.
In order to take the convolution of mono-energetic
coefficients, D*, with a Maxwellian energy distribution it is1041021010References[1] K. Ida et al., Phys. Rev. Lett. 91, 85003 (2003).
[2] S.P. Hirshman et al., Phys. Fluids 29, 2951 (1986).
[3] W.I. Van Rij and S.P. Hirshman, Phys. Fluids B 1, 563(1989).
[4] H. Maassberg et al., Phys. Fluids B 5, 3627 (1993).
[5] H. Maassberg et al., Phys. Fluids B 5, 3728 (1993).
[6] Y. Ogawa et al., Nucl. Fusion 32, 119 (1992).
[7] R. Kanno et al., Nucl. Fusion 37, 1463 (1997).
[8] A. Wakasa, S. Murakami et al., J. Plasma Fusion Res. 4,408 (2001).
[9] W. Lotz et al., Phys. Fluids 31, 2984 (1988)[10] C.M. Bishop, Rev. Sci. Instrum. 65, 1803 (1994)[11] S. Murakami et al., Nucl. Fusion 39, 1165 (1999).
[12] W.H. Press et al., Numerical recipes in C:The Art ofScinetific Computing, 2nd ed.

0-2-410 -6-4-20210 10 10 10 10ν*Fig. 4The normalized diffusion coefficient D * vs thenormalized collision frequency ν* at ρ = 0.5. The squaresymbols represent the results by DCOM with thetraining of the NNW and the solid lines represent theoutputs from NNW. The triangle symbols represent theDCOM results without the training of NNW and thedashed lines represent the results with NNW.

a) D ∗b) D∗

1053
101
1010-1
10-3
10-5

1053
101
1010-1
10-3-510

10-6-41010-2

ν∗Fig. 5100

102

10-5 10-4 10-2-3 10

10-110-610-410-2

ν∗G

100

102

0.2

0.4

0.6

0.8

1

ρ

The normalized monoenergetic coefficients, D* as a function D*(ν*, ρ, G). (a) ρ is constant (ρ = 0.5) and (b) G is constant (G =0.0).

4206