Application of an improved control variate technique to a local neoclassicaltransport code based on the δf Monte Carlo method
Seikichi Mattsuoka1 and Shinsuke Satake1
1National Institute for Fusion Science, JapanMay 15, 2012
1Introduction

Transport phenomena in a magnetically confinedplasma arises from two diﬀerent mechanisms: particle collision and turbulence. The former which is usually called
neoclassical (NC) transport determines an irreducible minimum transport level in a plasma. It also has a close relationship to the radial electric field, the bootstrap current, and so on which aﬀects the transport properties in aplasma through various mechanisms. Hence, the accurate
evaluation of the NC transport is a key issue to investigate the confinement performance of a device although theturbulence-driven transport often surpasses the NC transport.
The drift kinetic equation which involves the guiding
center motion of particles is a basis to investigate the NC
transport [1] and the particle simulation is widely used tonumerically solve it. Although it has a lot of advantages
such as its easiness to parallelize, applicability to multidimensional problems, etc., there is an inherentdrawback√
of a numerical noise proportional to 1/ N , where N denotes the number of markers used in the simulation. Theso-called δf Monte Carlo approach only solves the deviation of the total distribution function from the timeindependent part (usually Maxwellian) to reduce the noise.
This δf method can be interpreted as a control variatetechnique for the variance reduction in the Monte Carloframework [2].
We have developed a numerical transport simulation
code, FORTEC-3D [3, 4], based on the two-weight δfMonte Carlo method [5, 6]. FORTEC-3D have an advantage of including the finite orbit width (FOW) eﬀect in thedrift kinetic equation, which has been usually neglected inconventional numerical NC transport studies but becomes
important in high temperature plasmas due to the large
radial drift of trapped particles and the low collisionality.
It should be emphasized that the FOW eﬀect is expected
to become important in fusion reactors in the future. The
FOW eﬀect arises from the finite excursion of the guiding center orbits across the magnetic field lines and brings
higher order correction to the NC transport. The FOW effect on the NC transport and the formation of the Er has
been confirmed by comparing results of FORTEC-3D andthose of other codes which does not include it [4, 7]. The
codes used for the comparison above, however, are basedon the diﬀerent approaches to calculate the NC transport,thus more detailed and direct investigation for the FOW
eﬀect is required to quantitatively and individually assess
their influence on the NC transport. A new local NC transport code based on the same physical model as FORTEC3D except for the FOW eﬀect needs to be developed.

On the other hand, since such a local code isless time/CPU-consuming as compared to the original
FORTEC-3D due to its absence of the radial motion, it
can provide us to test a new improved control variate (CV)technique proposed by Kleiber, et al. in Ref. [8]. Their improved control variate method has been proved to reduce
the numerical noise drastically when applied to a simple
one-dimensional collisional model. Thus, to test the eﬀectof the new scheme at first, the new local code, which solves
the drift kinetic equation locally in four-dimensional phasespace, has been developed.
In this paper, we present initial results obtained by
our new local neoclassical transport code with the CVmethod. The rest of the paper is organized as follows:First, the two-weight δf Monte Carlo method with thenew CV method is briefly reviewed in Sec. 2. Section
3 devotes to numerical results with and without the newcontrol variate method. Finally, a summary is given inSec. 4.
2

δf Monte Carlo method

FORTEC-3D solves the drift kinetic equation for thedeviation of δf (z, t) = f − fM :dδf
∂t

≡
=[]()
∂+ v∥ + vd · ∇ − CTP δf∂t−vd · ∇fM ,(1)

where z denotes arbitrary phase space variables, f and fM
are the total distribution function and the Maxwellian fM ,v∥ = v∥ b and vd denote the parallel and drift velocity of
the guiding center particle, respectively, b is the normalized vector parallel to the magnetic field,and CTP is thepitch angle collision operator. For simplicity the new local
code at present only involves the pitch angle scattering andneglects the energy scattering for markers. In the Monte
Carlo method, it is required to solve the eq. (1) along thetrajectories of marker particles defined by ż, where dotdenotes the total time derivative, where the phase spacecoordinates, (ψ, θ, ζ, v∥ , v⊥ ), are used for the orbit calculation, and ψ, θ, ζ represent the radial label, poloidal√ andtoroidal angle variables respectively, and v⊥ =v 2 − v∥2and v is the magnitude of the velocity of the marker. Itis noted that ψ can be assumed constant in the local approach here due to the absence of the radial drift. For thetwo-weight δf scheme in FORTEC-3D, each marker has
two weights of w1 and w2 and the discretized marker particle distribution function is given in the extended phasespace F = F (z, w1 , w2 ; t). Then the solution of eq. (1),δf and the Maxwellian fM at any time can be evaluated0.003

obtained by the improved control variate method does not
reduce the noise compared to that obtained by the ordinary δf Monte Carlo method. This is accounted for by thesmall variance of w2 in the simulation. Since the noise inthe δf method mainly arises from increase in the varianceof w2 and a resultant weight-spreading in w1 , the noisein the ordinary δf method consequently does not becomelarge if the variance of w2 is does not increase. In fact,
the αi factor for the improved method is nearly kept unity
at all over the ξ-space, and α ≃ 1 represent the situationthat the ordinary δf method works well suﬃciently.

w/ CV
ordinary delta-f

0.002

δ f [a.u.]

0.001

0-0.001-0.002

4Summary and discussion-0.003-1-0.8-0.6-0.4-0.2

00.2

0.4

0.6

0.8

1

pitch angle ξFig. 1. The distribution function as a function of the pitch
angle variable, ξ. A red solid line represents a result obtained by using the improved control variate, and a blue
dashed line by the ordinary δf Monte Carlo approach.

∫∫using F as: δF = w1 F dw1 dw2 and fM = w2 F dw1 dw2 .
The time evolutions of w1,p and p2,p are given as follows:ẇ1,p=ẇ2,p=w2,p ∂fMψ̇,fM∂ψw2,p ∂fM
ψ̇,fM∂ψ−(2)(3)

where the p subscript denotes the index of the maker particle. It is noted that the v̇ = 0 is used because the potential
energy is assumed to be constant on a flux surface.
It is pointed out that the two-weight δf Monte Carlo
method described above does not fully exploit the advantage of the δf procedures [8]. Further reduction of
the noise is realized by introducing a new weight w̃1,p =cp − αw2,p , where cp = w1,p + w2,p and α is a numericalfactor specified later. With the new weight,
the total dis∫tribution f can be obtained by f = w̃1 F dw1 dw2 . The
parameter α is determined to minimize the variance of theresultant f by choosing
Cov [c, w2 ]α≡,V [w2 ](4)

where Cov[·, ·] and V[·] denote the covariance and the variance of given quantities. It is noted that α = 1 corresponds
to the ordinary δf Monte Carlo method and α = 0 to thefull-f method, respectively.
3

Numerical results

The new local code has been applied to a plasma of electron species with an axisymmetric (ζ-independent) magnetic field of which the magnitude is given by, B/B0 =1 − ϵt cos θ, where B0 correspond to the magnetic fieldstrength at the core and ϵt denotes the toroidicity. Here,B0 = 1.88T and ϵt = 0.05 are used. The number of marker
particles used is 32, 000.
To see how the improved control variate scheme reduces
the numerical noise, the δf part of the distribution is investigated. Figure 1 shows resultant δf at the steady state
as a function of the pitch angle ξ = v∥ /v. It is noted that
α is evaluated only in the ξ direction, that is, the weights
w1 and w2 are integrated in the θ and v directions to evaluate the eq. (4). As can be seen in this figure, the result

In this paper, we report the application of an improved
control variate technique for a local neoclassical transport
code, which aims a direct comparison with our FORTEC3D code to investigate the FOW eﬀect in more detail.
The eﬀect of the variance reduction has been explored by
the introduction of the improved control variate to theordinary two-weight δf Monte Carlo method. The improved control variate is calculated only in the ξ-space
although the drift kinetic equation solved in the code isfour-dimensional. The result has shown that no significant improvement has been observed when compared tothe ordinary δf method. This has occurred from the factthat the ordinary δf method has already worked well dueto the small variance of the w2 .
The variance of w2 in general increases as the radialdrift of particles increases since the evolution equation is
proportional to the ψ̇ which represent the drift velocity.
Thus, the control variate is expected to be more eﬀective if
we apply it to a higher temperature and/or less collisional
plasma, and we will report in the presentation the resultsof parameter survey calculations for the improved control
variate. In addition, the multi-dimensional control variate
method and the energy scattering for the marker particles
are also implemented and results will be presented.
Acknowledgment
This work is supported in part by JSPS Grantin-Aid for Young Scientists (B), No.
23760810 andin part by the NIFS Collaborative Research Program,
NIFS12KNST038.
References[1] P. Helander and D. J. Sigmar., Collisional Transport in
Magnetized Plasmas, Cambridge University Press, Cambridge, 2002.
[2] A. Y. Aydemir., Physics of Plasmas, 1, 822, 1994.
[3] Shinsuke Satake, et al., Plasma and Fusion Research, 3,S1062, 2008.
[4] Seikichi Matsuoka, et al., Physics of Plasmas, 18, 032511,2011.
[5] S. Brunner, E. Valeo, and J. A. Krommes. Physics of Plasmas, 6, 4504, 1999.
[6] W. X. Wang, et al.,Plasma Physics and Controlled Fusion,
41, 1091, 1999.
[7] Seikichi Matsuoka, et al., Plasma and Fusion Research, 6,1203016, 2011.
[8] R. Kleiber, et al.,Computer Physics Communications, 182,1005, 2011.