Proceedings of ITC/ISHW2007

Computations of neoclassical transport in stellarators usinga δ f method with reduced variance ∗Klaus ALLMAIER1) , Sergei V. KASILOV1,2) , Winfried KERNBICHLER1) , Georg O. LEITOLD1) ,Viktor V. NEMOV1,2)
1)

Association EURATOM-ÖAW, Institut für Theoretische Physik - Computational Physics,TU Graz, Petersgasse 16, A–8010 Graz, Austria
2)Institute of Plasma Physics, National Science Center “Kharkov Institute of Physics and Technology”,
Akademicheskaya Str. 1, 61108 Kharkov, Ukraine
An improved δ f Monte Carlo method for the computation of neoclassical transport coefficients in stellarators
is presented. Compared to the standard δ f method without filtering, the computing time needed for the same
statistical error decreases by a factor proportional to the mean free path to the power 3/2.
Keywords: bootstrap current, stellarator, δ f Monte Carlo, variance reduction
A standard δ f Monte Carlo method for the computation of neoclassical transport coefficients [1] assumes solution of the linearized drift-kinetic equation taking into
account the source term in the equation for test particle weights which evolve with time. This method has agood convergence for tokamaks where the variance in alltransport coefficients including bootstrap coefficient has nostrong dependence on plasma collisionality. However, instellarators, the variance in the bootstrap coefficient increases for this method as a square of the mean free pathdue to the accumulation of large random contributions to
the test particle weights which occur in the phase space region occupied by trapped particles. In Ref. [2] a method
has been presented which is a combination of the standard
Monte Carlo method with a method employing constantparticle weights and re-discretizations of the test particle
distribution in phase space. There, this method has beentested in confinement regimes with negligible radial electric field. Below this method is described in more detail for
general confinement regimes.
Mono-energetic transport coefficients are determinedby the steady state solution of the linearized drift kinetic
equation for the normalized perturbation of the distributionfunction fˆ (marker)!
∂LD fˆ ≡+ Vg · ∇ − LC fˆ = ψ̇ ≡ Vg · ∇ψ, (1)∂twhere LC , Vg , and ψ̇ are Lorentz collision operator, drift
velocity and its co-variant ψ-component, respectively, andvelocity space variables are total energy and perpendicular∗ This work, supported by the European Communities under the contract of Association between EURATOM and the Austrian Academy ofSciences, was carried out within the framework of the European Fusion
Development Agreement. The views and opinions expressed herein do
not necessarily reflect those of the European Commission. Additional
funding is provided by the Austrian Science Foundation, FWF, under contract number P16797-N08.
author’s e-mail: allmaier@itp.tugraz.atadiabatic invariant. Here, ψ is a flux surface label, and amarker is defined through the local Maxwellian distribution function f M and the total distribution function f viaf = f M − fˆ∂ f M /∂ψ. The mono-energetic radial diffusion coefficient and the normalized bootstrap coefficient,
respectively, are given byDmono=−* Z1+1ˆψ̇ ,dλfh|∇ψ|i2 21(2)

−1

λbb+
* Z131ˆ= −dλ f λB ,ρL B0 h|∇ψ|i 2(3)−1where λ = vk /v is the pitch parameter, ρL is the Larmor
radius in the reference magnetic field B0 , B is the magnetic field module and h. . .i denotes the average over thevolume between neighboring flux surfaces. The quantity λbb is linked to the equilibrium (bootstrap) current
density jk and gradient of the pressure p by λbb =− jk B (c h|∇ψ|i dp/dψ)−1 . In the following Dmono is normalized√ by the plateau diffusion coefficient Dplateau =πvρ2L (8 2ιR)−1 where ι and R are the rotational transformand major radius, respectively.
In order to introduce the Monte Carlo operator it isconvenient to re-write (1) in the integral form using aGreen’s function G defined byLD G(t, z, z0 )=G(0, z, z0 )=(4)0(g(z0 ))−1/2

δ(z − z0 ),(5)where z = (ϑ, ϕ, λ) and g is a metric determinant of flux
coordinates (ψ, ϑ, ϕ). This Green’s function is normalized
to 1,Zd3 z (g(z))1/2 G(t, z, z0 ) = 1.

(6)

Proceedings of ITC/ISHW2007
Thus, a formal solution to Eq. (1) isZ3
1/2 ˆf (t, z) =d z0 (g(z0 )) G(t − t0 , z, z0 ) fˆ(t0 , z0 )Zt×dt0G(t − t0 , z, z0 )ψ̇(z0 )

(7)

t0

If a steady state solution is looked for, fˆ(t, z) =ˆf (z), Eq. (7) becomes an integral equation for F(z) =(g(z))1/2 fˆ(z) given below also in the operator form,ZF(z) =d3 z0 K(z, z0 )F(z0 ) + Q(z) ≡ K F + Q, (8)where K(z, z0 ) = (g(z))
time step and

1/2∞
Xw0 λk B(zk )=k=0G(∆t, z, z0 ), ∆t is the integrationZ∆t31/2

d z0 (g(z)g(z0 ))dt0G(t0 , z, z0 )ψ̇(z0 )=Wk0≈ (g(z))1/2 ψ̇(z)∆t.

∞
X=K k Q = C0(9)

k=0

zk=∞
Xw0 δ (z − zk ),

(10)

k=0

w0 = ψ̇(z0 )∆t,Z(∆t, zk−1 ),where C0 = d3 z (g(z))1/2 and the random starting point
z0 is chosen with the probability density δ (z − z0 ) =C0−1 (g(z))1/2 . The averages (2) and (3) are given by expectation values asRDmono=−λbb=−1∞Xh|∇ψ|i2 k=0w0 ψ̇(zk ),

(11)

∞X
3w0 λk B(zk ).
ρL B0 h|∇ψ|i k=0(13)1 X
Wk λk B(zk ),K→∞ K
k=0limkX

w j.

(14)

(15)

j=0

Introducing the Monte Carlo operator, Z(∆t, z0 ), being therandom position of a test particle starting at z0 after a single time step modeled in a standard way [4] by the randomchange of λ in accordance with LC and integration of particle drift equations, the kernel of the integral equation isgiven by an expectation value K(z, z0 ) = δ (z − Z(∆t, z0 )).
The solution of (8) by direct iterations can be presented as an expectation value of an integral along thestochastic orbit,
Flim Wk λk B(zk )k→∞

K=Z=Q(z)are the connection length and mean free path, respectively,
contribute to λbb . In addition, the contribution of a particle from the boundary layer is ∆λ times smaller than of anormal passing particle because of a higher trapping probability. Therefore, the variance in λbb scales for this methodas (lc /Lc )2 in the long mean free path regime.
It should be noted that distribution of test particles after each step remains to be the equilibrium distribution,
δ (z − zk ) = C0−1 (g(z))1/2 .
Therefore w0 λk− j B(zk− j ) = w j λk B(zk ) where w j =ψ̇(z j )∆t, and(12)When k∆t exceeds a few collision times, the correlation
between zk and w0 is lost and, therefore, such terms in (11)tend to zero, e.g. w0 ψ̇(zk ) → Rw0 ψ̇(zk ) = 0 because theexpectation value w0 = C0−1 ∆t d3 z (g(z))1/2 ψ̇(z) = 0 due
to Liouville’s theorem. Thus, a finite sum over k is sufficient in (11). The method of constant test particle weights
described by (11) and (12) has rather low variance in computations of Dmono , however, for λbb variance has a very
unfavorable scaling with collisionality. Indeed, only theorbits originating in the boundary layer in velocity space of
the width ∆λ ∼ (Lc /lc )1/2 where Lc = 2πR/ι and lc = vτc

The procedure in (13) corresponds to a standard δ fmethod [1] where the test particle weight Wk is an integralof ψ̇ along a stochastic orbit. In a tokamak, the variance in
λbb does not scale with the collisionality, and the required
CPU time for this method scales linearly with lc /Lc . However, in stellarators variance in λbb again recovers the scaling (lc /Lc )2 because due to the non-zero bounce-averaged
drift of trapped particles large contributions to Wk are acquired, which scale as ψ̇τc and which become weakly correlated with the values of λk after detrapping of test particles. To overcome this problem, in Ref. [5] trapped particles with large Wk are replaced with particles with Wk = 0(large weights are filtered out) which formally introducessome bias in the result.
For a formally “unbiased” method it is convenient tosplit the source in (8) into “passing” and “trapped” sources
Q p = χQ and Qt = Q − Q p using
1
1 + tanh (λ − λt−p )/∆λ ,(16)2where λt−p is a trapped-passing boundary, and solve theproblem with each source independently. Results for transport coefficients for these two sources are added up atthe end. The problem with Q p is solved with the standard method (13) because accumulation of large weights isavoided there. For the problem with Qt the formal solutionto (8) is presented as F = F M + ∆F M where F M satisfies
an equation which differs from (8) only by a source term,
χ=F M = K F M + QM ,

(17)

where
M−1

QM

=

∆F M

=1 X k
K Q,
M k=0
M−1
X
k=0

1−!
k+1
K k Q.
M(18)(19)

Proceedings of ITC/ISHW2007
In order to derive (17) the original equation (8) is presentedasm−1
X
F =K k Q + K m F,(20)10log(N)8k=0m=0

K k QM ,42

m=1Denoting the last term in (21) with F M and substituting
there F in the form of the series (10) one obtains
 M−1MX
1  X mm
K Q+K F FM ≡(22)
M m=0
m=1∞
X

[a.u.]

6

where m is an arbitrary natural number. Averaging the r.h.s.
of (20) over 1 ≤ m ≤ M yields
 M−1
MX

1  X m

K Q+(21)K m F F = ∆F M +M=λbb0020406080100
number of iterations [a.u.]

120140Fig. 1 Bootstrap coefficient λbb (green) and logarithm of thenumber of simulation particles (red) plotted over thenumber of iterations.

(23)

k=0

610tCPU ⋅ σ2 [a.u.]which is a series solution to Eq. (17).
Equation (17) describes one iteration of the solution
procedure. Given the source term Q such iteration provides
Q M which is used as a source term Q for the next iteration.
Within an iteration, each of the test particles performs M−1steps, and the quantity Q M is computed by scoring weights
on the 3D grid (except the first iteration). The contribution of the iteration to the distribution function, ∆F M , isformally used in the averages (2) and (3) so that they arecomputed directly in analogy to (11) and (12). For the first
iteration with Q being the original source (9) the algorithmwith alternating weights (15) is used. After the first iteration these test particle weights are divided by M so that test
particles represent the discretized distribution Q M used asa source within the second iteration. Starting from this second iteration, an algorithm with fixed test particle weights(and scoring Q M on the grid) is used. After the second iteration, the module of test particle weight is fixed. Due toannihilation of the weights on the grid and fixed moduleof the weight, the number of test particles needed for sampling the source term from the grid is decreasing with iterations and iterations are stopped when this number is below
one, (see Fig. 1). The number of steps M for a single iteration is chosen to be much smaller than collision time andlarge enough in order to fill the grid using a limited numberof test particles. Since source terms generated in this way
are small in the passing and boundary region, particles aregenerated there with smaller weights and particles which
enter the boundary layer from the trapped side are split insuch a way that the number of test particles in the passing
and trapped regions is of the same order. As a result, variance of this method is reduced to the scaling lc /Lc . In addition, due to the decay of test particle number with iterations, the CPU time is also reduced and scales as (lc /Lc )3/2
for a given accuracy (see Fig. 2) which is much better than
the scaling (lc /Lc )3 of the standard δ f method without a filter. Results of testing of the method for a few toroidal devices stay in good agreement with results of NEO-2, a field

8104102

10010−510−410−3

Lc/ lc10−210Fig. 2 CPU-time multiplied with the variance of the bootstrap
coefficient σ2 plotted over the collisionality parameter(red). The green line shows the scaling (lc /Lc )3/2 .

line tracing code which computes transport coefficients inarbitrary collisionality regimes [3], as shown in Figs. 3-6.
Benchmarking of the results from computations with finite
radial electric fields with other methods has been started,first results can be found in Ref. [6].
[1] M. Sasinowski and A.H. Boozer, Phys. Plasmas 2, 610(1995).
[2] K. Allmaier, S. V. Kasilov, W. Kernbichler, G. O. Leitold
and V. V. Nemov, 34th EPS Conference on Plasma Physics,Warsaw, 02–06 July 2007, P-4.041, (2007).
[3] W. Kernbichler, S. V. Kasilov, G. O. Leitold, V. V. Nemov,
and K. Allmaier, 33rd EPS Conference on Plasma Physics,Rome, 19–23 June 2006, ECA Vol. 30I P-2.189, (2006).
[4] A. H. Boozer and G. Kuo-Petravic, Phys. Fluids 24, 851(1981).
[5] M. Yu. Isaev, S. Brunner, W. A. Cooper, et al., Fusion Science & Technology 50, 440 (2006).
[6] K. Allmaier, C.D. Beidler, M.Yu. Isaev, S.V. Kasilov, W.
Kernbichler, H. Maassberg, S. Murakami, D.A. Spong, V.
Tribaldos, ICNTS - Benchmarking of Bootstrap Current
Coefficients, this workshop, (2007).

Proceedings of ITC/ISHW2007

1210110Dmono/Dplateau

Dmono/Dplateau

10010010−110−110−410−3−21010Lc/ lc−110010−410Fig. 3 Normalized diffusion coefficient Dmono /Dplateau for LHD
with R=375cm vs. collisionality parameter Lc /lc at half
plasma radius computed by NEO-MC (blue) and NEO-2(red) for Er /v/B = 0 (circles), 3 · 10−5 (stars), 1 · 10−4(squares), 3 · 10−4 (diamonds), 1 · 10−3 (triangles)−310−210Lc/ lc−110010Fig. 5 Normalized diffusion coefficient Dmono /Dplateau for W7-X
standard configuration vs. collisionality parameter Lc /lcat half plasma radius. Markers and colors are the same asin Fig. 3121100.8

80.6
λbb

λbb.

640.4
0.2

200−0.2−410−310−210Lc/ lc−110010Fig. 4 Normalized bootstrap coefficient λbb for LHD with
R=375cm vs. collisionality parameter Lc /lc at halfplasma radius. Markers and colors are the same as inFig. 3−410−310−210Lc/ lc−110010Fig. 6 Normalized bootstrap coefficient λbb for W7-X standard
configuration vs. collisionality parameter Lc /lc at half
plasma radius. Markers and colors are the same as inFig. 3