Gyrokinetic Simulations Using Particles

1
Gyrokinetic Simulations Using Particles

• L. Villard1
• with many thanks to
• S. Brunner1, P. Angelino2, A. Bottino3, R.
  Hatzky4, Y. Idomura5, S. Jolliet5,
• B.F. McMillan1, T.M. Tran1, T. Vernay1
• (1) Centre de Recherches en Physique des Plasmas,
  Association Euratom - Suisse, EPFL, 1015
  Lausanne, Switzerland
• (2) Association Euratom - CEA/DSM/IRFM,
  Cadarache, France
• (3) Max-Planck Gesellschaft, Association Euratom
  - IPP, Garching, Germany
• (4) Computer Center, IPP, Max-Planck
  Gesellschaft, Garching, Germany
• (5) Japan Atomic Energy Agency, Taitou, Tokyo
  110-0015, Japan

2
ITER the way to fusion
SourceITER

• EU, Japan, USA, China, India, South Korea,
  Russian Federation
• gt 5.3G construction cost (1/3 G /y 0.06 of
  world overall RD)
• Virtually inexhaustible, environmentally benign
  source of energy
• from Deuterium and Lithium gives Helium.
  (Tritium is recycled)

3
Magnetic confinement tokamak
Trapped particle
Larmor radius rL
Passing particle
particle trajectory
4
Complexity many nonlinearities
Neutral Beams
heat

• Geometry of magnetic configuration is an
  essential feature of fusion plasma physics

5
Timescales in the ITER plasma
machine lifetime
energy confinement
ion cyclotron
electron cyclotron
1 shot
turbulence

• Physics spans several orders of magnitude
• Direct Numerical Simulation (DNS) of everything
  is unthinkable
• Need to separate timescales using approximations
• How to integrate all phenomena in a consistent
  manner?

6
Kinetic effects wave-particle interaction

• General distribution function in 6D phase space
• To be solved with consistent electromagnetic
  fields

7
Outline

• 1. Turbulence and transport, GyroKinetic theory
• 2. Numerical schemes to solve GK equations
• 3. Critical issue no.1 numerical noise
• 4. Critical issue no.2 geometry, anisotropy
• 5. Critical issue no.3 long, (quasi-)steady
  state simulations
• 6. Collisions in PIC
• 7. High Performance Computing issues
• 8. Outlook and conclusions

8
Outline

• 1. Turbulence and transport, GyroKinetic theory
• gyrokinetic modelling of fusion plasmas
• 2. Numerical schemes to solve GK equations
• 3. Critical issue no.1 numerical noise
• 4. Critical issue no.2 geometry, anisotropy
• 5. Critical issue no.3 long, (quasi-)steady
  state simulations
• 6. Collisions in PIC
• 7. High Performance Computing issues
• 8. Outlook and conclusions

9
1. Turbulence and Transport

• Hot plasma core, colder edge gt grad T
• grad T gt (grad T)crit gt instabilities, even when
  macroscopic (MHD) modes are stable
• small scales, i.e. spatial scales rs (Larmor
  radius)
• large scales, i.e. spatial scales system size a
• low frequency, i.e. w ltlt Wi (cyclotron frequency)
• Growth, then nonlinear saturation of modes
• gt turbulent state
• gt transport (heat, particles, momentum)
• anomalous, essentially by collisionless
  processes
• Turbulent transport gtgt collisional transport
• In spite of this, typical grad T 106K / cm in
  tokamaks
• Need a theory and numerical tools!

10
Boltzmann/Vlasov Maxwell kinetic model

• Boltzmann equation for single particle
  distribution function fs(q,p)

Along characteristics (orbits), determined by
particle equations of motion
.,. Poisson bracket in canonical coordinates
(q,p)
C(fs , fs) collision operator
Hs(q, p) Hamiltonian of collisionless single
particle motion
Moments of the distribution function
Maxwells equations
11
Spatio-temporal scales in a magnetic fusion
plasma and range of validity of Boltzmann/Vlasov,
Gyrokinetic and MHD models
Y. Idomura et al., C.R. Physique 7 (2006) 650-669
12
Gyrokinetic model basic assumptions

• Small parameter eg , with

characteristic length scale of equilibrium B

• Small parameter

13
Gyrokinetic model

• Assume
• Average out the fast motion of the particle
  around the guiding center
• Fast parallel motion, slow perpendicular motion
  (drifts)
• Strong anisotropy of turbulent perturbations (//
  vs perp to B)

phase space dimension reduction 6D ---gt 5D
14
Guiding center motion in a tokamak
Passing (top) and trapped (bottom) particle
orbits the guiding centres almost follow the
equilibrium magnetic field. Left 3D view.
Right projection on the poloidal plane
Y. Idomura et al., C.R. Physique 7 (2006) 650-669
15
Guiding center coordinates and Gyro-centre
transform
Guiding center coordinates
generalised parallel velocity
magnetic moment
16
Gyrokinetic equations
It is an advection equation along nonlinear
characteristics
curvature and grad-B
magnetic and ExB nonlinearities
parallel motion (FAST)
drifts (SLOW)
parallel velocity nonlinearity
mirror term
conservation of the magnetic moment
generalized potential
17
Gyrokinetic perturbed field equations

• Poisson (or quasi-neutrality) equation, here
  with Boltzmann electrons, linearized ion
  polarization density, long wavelength approx.

gyro-center perturbed ion density
flux-surface-averaged potential

• Ampères law, here neglecting dB// and expanding
  

gyro-center perturbed ion and electron currents
Integral partial differential system of
equations, inhomogeneous, linear
18
Properties of GK-Maxwell equations (in the
collisionless, sourceless case)

• The distribution function f (and any arbitrary
  function of f ) is conserved along nonlinear
  characteristics
• The phase space volume is conserved
• High frequencies (e.g. Langmuir and cyclotron
  waves) are eliminated
• large time steps are allowed
• Conservations of particle number, momentum and
  energy
• can serve as a stringent test of the quality of
  numerical simulations

19
Gyrokinetic-Maxwell summary
Advection in incompressible phase space
Field equations Poisson, Ampere
If collisions and sources are included
df/dtC(f)S
20
Gyrokinetic-Maxwell summary

• A time evolution partial differential equation
  (PDE) describing the advection in 5D phase space
  of the distribution function f
• A set of nonlinear, coupled 5 ordinary
  differential equations (ODEs) for the
  characteristics
• A (usually linearized) set of integral partial
  differential equations for the perturbed fields
  (potentials f , A// ) in 3D real space

21
Outline

• 1. Turbulence and transport, GyroKinetic theory
• 2. Numerical schemes to solve GK equations
• Lagrange Particle in Cell (PIC)
• Particle (markers) Klimontovich and Finite
  Element representations
• Fourier filter
• Real space filter
• Other schemes Eulerian, Semi-Lagrangian
• 3. Critical issue no.1 numerical noise
• 4. Critical issue no.2 geometry, anisotropy
• 5. Critical issue no.3 long, (quasi-)steady
  state simulations
• 6. Collisions in PIC
• 7. High Performance Computing issues
• 8. Outlook and conclusions

22
2. Solving GK numerical methods
Solving the GyroKinetic - Maxwell system of
equations requires the following elements

• Time-advancing scheme for a distribution function
  in 5D phase space
• e.g. Runge-Kutta, Verlet, leapfrog,
• ODE integrator (for characteristics)
  (alternatively, a finite difference scheme in 5D
  phase space)
• Field solver for integral-differential operator
  in 3D real space
• e.g. Fast Fourier Transforms, finite differences,
  finite elements,

23
Solving GK numerical approaches

• Lagrangian Particle In Cell (PIC), Monte Carlo
• Sample phase space (markers) and follow their
  orbits
• Noise accumulation is difficult to control in
  long simulations
• Semi-Lagrangian
• Fixed grid in phase space, trace orbits back in
  time
• Multi-dimensional interpolation is difficult
• Eulerian
• Fixed grid in phase space, finite differences for
  operators
• CFL condition overshoot versus dissipation

• Common to all three Poisson (Ampere) field
  solver
• Various methods finite differences, finite
  elements, FFT,

Developing several approaches is not a luxury
cross -comparisons are a necessity.
24
Lagrange PIC essential steps
grid cell
25
Lagrange PIC basic approach

• Birdsall C.K. and Langdon A.B., Plasma physics
  via computer simulations, IOP Series in Plasma
  Physics (Taylor Francis, 2004)
• Simple, 1D electrostatic problem

sum over physical particles
26
Lagrange PIC basic approach

• Use N superparticles (N ltlt Nphysical)
• Define an Eulerian grid xj, j1..Nx on which
  fields are computed
• Define a temporal grid tk, k1..Nsteps

Essential step represent the charge on the
Eulerian grid from superparticle positions and
charges. This is called charge assignment
Contribution of one superparticle
27
Lagrange PIC basic approach

• In the previous example, first-order weighting
  has been applied

• Equivalent particle cloud picture

Hence the name particle-in-cell
28
Lagrange PIC basic approach

• Solve Poissons equation, e.g. with finite
  differences

• Obtain electric field at particle position using
  same weighting as for the charge assignment

• Higher order weightings can be applied

29
Lagrange-PIC finite element
particle representation
finite element representation
charge assignment
0-th and 1st order
quadratic
cubic
Examples of B-spline shapes
30
Lagrange-PIC finite element

• B-splines have interesting properties

• The sum of all spline elements at any point x is
  unity

Useful for conservation properties

• They have unitary surface

• The derivative of the spline shape of order n can
  be obtained interms of the spline shape of order
  n-1

Useful to obtain E from f

• Their Fourier transform in a periodic system of
  Nx intervals is

Useful for statistical noise reduction
31
Lagrange-PIC finite element

• Field equation variational Galerkin formulation
• Mutiply Poissons equation with test function
  h(x) and integrate over the spatial domain
• Integrate by parts
• Represent h(x) on the same basis functions as for
  f(x)
• Discretized variational problem
• Linear algebraic system

for all h(x)
for all h(x)
32
Lagrange-PIC finite element
With the particle representation
we obtain
What is the corresponding finite element
representation of r(x) ?
Let us call this representation r fe (x)
The variational Galerkin formulation applied to
the identity problem
gives
Mass matrix

• The mass matrix inversion results in some
  problems for the spline representation of the
  density, e.g. non-positiveness (see further)

33
One particle (one delta function) on F.E.
r fe d(x-xp)
r fe
r PIC M r fe

• The finite element representation introduces
  grid-related unphysical oscillations. Positivity
  is NOT ensured.
• The field equations are linear, so this will
  simply add up for each marker! Only in the limit
  Np infinity will these oscillations tend to zero.

34
Application of a Fourier filter

• Numerical statistical noise is slow to converge
  (1/Np1/2)
• Error of overshoots due to mass matrix inversion
  is large
• Error can be greatly reduced by the application
  of filters
• Fourier

Apply filter coefficients Fk

• Real space filter, e.g. mass matrix filter

35
Effects of Fourier filter
Illustration with representation of r sin (2p
mx)
m1...7
r (x)

• NB with only 2 markers/cell, filtered results
  are excellent!
• 36 markers / mode in the filter
• The statistical noise thus depends on the number
  of markers per Fourier mode retained in the
  filter, not on the number of markers per cell

36
Efficiency of Fourier filter
Compare Fourier-filtered Np/Nmodes36 and
unfiltered with Np/Nx36

• Superiority of Fourier-filtered even with 20
  times less markers!

37
Euler (Vlasov) (continuum)
38
Semi-Lagrange
39
ITER global GK turbulence simulation
requirements

• Normalized size system size / gyration radius
  1 / r a / rs rs cs/Wi cs(Te/mi)1/2
• ITER 1 / r 1000
• Simplest model electrostatic, GK ions, adiabatic
  electrons
• Ion Temperature Gradient (ITG) driven
  turbulence
• Relevant scales resolve up to krrs, kqrs lt 1
• 3D field solver min 4 pts/wavelength
• Nr (2/p) 1/r 640, Nq 4 /r 4000, Nj
  Nq/q 3000
• Phase space resolution 300 markers / cell ()
  
• or Nm8, Nv//32
• () Needs to be assessed by numerical convergence
  studies
• Time until statistical steady-state 2000 a /
  cs ()
• () Needs to be assessed by series of (long)
  runs with independent initial conditions

40
ITER global GK turbulence simulation
requirements

• Field solver 600 x 4000 x 3000 7.2 x109 grid
  points.
• CPU cost scaling (1/r)3
• 2000 x 109 particles (Lagrange-PIC).
• CPU cost scaling (1/r)3
• 2000 x 109 grid points (Eulerian /
  semi-Lagrange).
• CPU cost scaling (1/r)3
• Number of time steps, assuming Wi Dt 1, Nsteps
  t/Dt 2000 (a/cs) Wi 2000 (1/r)
  2x106
• CPU cost scaling (1/r)
• The computational cost scales as (1/r)4 and is
  prohibitive for ITER-size parameters
• Improvements must be found !

41
Outline

• 1. Turbulence and transport, GyroKinetic theory
• 2. Numerical schemes to solve GK equations
• 3. Critical issue no.1 numerical noise
• Statistical interpretation of PIC
• Control Variates
• Loading issues
• 4. Critical issue no.2 geometry, anisotropy
• 5. Critical issue no.3 long, (quasi-)steady
  state simulations
• 6. Collisions in PIC
• 7. High Performance Computing issues
• 8. Outlook and conclusions

42
3. Critical issue no.1 numerical noise
Lagrange Particle-In-Cell (PIC)
err
field solver
get electric field

• Problem sampling error accumulation. Possible
  cures
• a) control variates
• b) importance sampling
• c) physics based filter
• d) explicit noise-control source operators

43
3.a) Statistical aspect of Lagrange-PIC

• Obtaining the density and currents from the
  markers, i.e. obtaining moments of the
  distribution funciton, can be interpreted as a
  Monte Carlo integration

• Even small errors in the evaluation of these
  moments can cause a systematic corruption of the
  simulation in a realtively short period of time.
  S.J. Allfrey, R. Hatzky, CPC 154 (2003) 98

See also A.Y. Aydemir, Phys. Plasmas 1 (1994) 822
44
3.a)

• Let p(Z) the probability density function (PDF)
  of N marker positions in phase space

Expectation value of gn
Error ebv in this estimate
where sn2 is the variance in this estimating
function gn
with
importance sampling choose
Nonzero variance due to finite spatial
localisation of Ln
45
3.b) Control variates df method
Only df is sampled with markers
marker weight
0 if f0 is a fct of constants of unperturbed
motion
We obtain the time evolution equation for the
weights
46
3.b) control variates df method

• Improve estimate for bn, i.e. reduce the variance
  sn2

df
The error in this estimate is proportional to the
sqrt of the variance of the function
instead of
Thus the error is reduced by a factor of the
order of
47
3.b) df method (continued)
Only the perturbed part of Poisson
(quasi-neutrality) equation is solved

• Sampling noise is reduced if

• The smaller the fraction of f is sampled with
  markers, the lower the sampling noise will be
• The basic idea is to replace as much as possible,
  in the field equations, Monte Carlo (markers)
  integration with analytical expressions
• Enhanced control variates

Effective if the non-Bolzmann part is much
smaller than the Boltmann part
48
Loading issues

• Idea optimize the pdf of initial loading
  positions Zp(t0) so as to mimimize the variance
  of the weights at later times
• Run a 1st simulation until ttload
• Obtain the pdf of weight distribution, traced
  back to the Zp(t0) positions
• Run a 2nd simulation using this pdf for the new
  loading

Physics based filter

• Supress obviously unphysical modesFourier
  B-field aligned filter see later

Noise control source operators

• Aim at long, quasi-steady, statistically
  converged simulations see later

49
3.c) Numerical convergence loading issues

• Which random number generator to use to load
  the initial marker positions?
• Pseudo random generators
• Linear congruential, Twister-Mersenne,
• give a convergence 1/Np1/2 for Monte Carlo
  integration (CLT applied to variance of sum)
• For 5D case, this is competitive against
  grid-based methods that have a convergence worse
  than h5/2
• Can we do better? Yes, with quasi-random
  generators!

50
Loading using quasi-random generators

• Haltons sequence
• Let r be a prime number
• Let i be a natural number written in base r
• The i-th number of Haltons sequence is
• Gives close to uniform distribution in 0,1

51
Pseudo-Random vs Quasi-Random

• Haltons in different bases for different
  dimensions give better coverage of phase space

52
Pseudo-Random vs Quasi-Random
N104 32 bins

• Quasi-random (Haltons) give a better restitution
  of the desired (here uniform) probability
  distribution function.

53
Pseudo-Random vs Quasi-Random

• Quasi-random (Haltons) give better accuracy in
  integration and higher order of convergence in
  error, i.e. almost 1/Np

54
ITG gyrokinetic global linear simulation

• CYCLONE base case A.M. Dimits et al, PoP 7
  (2000) 969
• Parameters inspired by experiment on DIII-D
  tokamak
• 1/r180

55
ITG gyrokinetic global linear simulation
Contours of perturbed potential
Linearised df

• CYCLONE base case 1/r180. GYGLES code M.
  Fivaz et al, CPC 111 (1998) 27

56
ITG gyrokinetic global linear simulation
Convergence of the growth rate with number of
markers Np
Halton sequence loading

• Pseudo-Random loading convergence 1/Np1/2
  (left)
• Halton sequence loading convergence 1/Np
  (right)

57
ITG gyrokinetic global linear simulation

• ITER-size cyclone case 1/r 1120

• m170 wavelengths, solved with Nq64 grid points
  how is it possible?

58
ITG gyrokinetic global linear simulation

• Checking convergence with Np
• ITER-size cyclone case 1/r 1120

Growth rate
Std deviation of growth rate
59
ITG gyrokinetic global linear simulation

• Checking the power balance

(1/2E) dE/dt
(1/2E) j.E
Instantaneous growth rate vs time
60
Outline

• 1. Turbulence and transport, GyroKinetic theory
• 2. Numerical schemes to solve GK equations
• 3. Critical issue no.1 numerical noise
• 4. Critical issue no.2 geometry, anisotropy
• Straight-field-line coordinates
• Field-aligned Fourier filter
• Benefits for noise problem
• Benefits for field solver problem
• Beware of simplified geometrical models
• 5. Critical issue no.3 long, (quasi-)steady
  state simulations
• 6. Collisions in PIC
• 7. High Performance Computing issues
• 8. Outlook and conclusions

61
4. Critical issue no.2 geometry, anisotropy
Small scales perpendicular to B
Next page
Color contours of perturbed electrostatic
potential
Long scales parallel to B

• Gyrokinetic ordering
• High k//rs modes are unphysical gt filter them
  out!

62
Straight-field-line coordinates
Poloidal angle-like coordinate q
toroidal angle j

• Choose q such that B field lines are straight in
  (q,j) plane
• In these coordinates the expression of the
  parallel wave number is simpler

63
Field-aligned Fourier Filter (FAFF)

• From gyrokinetic ordering, k//rs r ltlt 1
• Short wavelengths along B are unphysicalgt filter
  them out!

Writing the perturbations in Fourier space
The physical perturbations are thus such that
The number of poloidal Fourier modes to be
retained is thus 10, independently of the system
size !
64
Field-aligned Fourier Filter (FAFF)

• With FAFF, the scalability of computational cost
  with system size is dramatically improved
• Constraint on timestep from (fast) parallel
  motion

Without FAFF, Dmma/rs1/r, leading to
With FAFF, DmqR/a, leading to
With FAFF, reduction of computational cost by a
factor
For ITER-size, this is a factor 200 !
65
Field -aligned Fourier filter (FAFF)
1/r40, R/LT12, R/a5, Np16x106, grid
128x128x64

• Turbulent perturbations. Detail on a magnetic
  surface.

S. Jolliet, et al., CPC 2007
Without FAFF unphysical small scales in the
parallel direction develop
With FAFF unphysical small scales in the
parallel direction are suppressed
66
Fighting noise with field -aligned Fourier filter
(FAFF)

• Heat flux versus time

Without FAFF
S. Jolliet, et al., CPC 2007
With FAFF
Without FAFF, the noise increases the heat flux
67
Fighting noise with field -aligned Fourier filter
(FAFF)

• Checking energy conservation

Without FAFF
S. Jolliet, et al., CPC 2007
With FAFF Dm7
Without FAFF, energy non-conservation rises
indefinitely
68
Fighting noise with field -aligned Fourier filter
(FAFF)

• Energy of toroidal modes

Without FAFF
With FAFF, Dm7
S. Jolliet, et al., CPC 2007
Without FAFF, there is an unphysical growth of
ITG perturbations and of Zonal Flows
69
Markers and finite element representations
The marker-sampled perturbed density is
Field quantities are represented on finite
elements
Simplest case, 1D
Galerkin formulation
Mass matrix M
b
70
ITER global GK turbulence simulation
requirements, with field-aligned Fourier

• Normalized size system size / gyration radius
  1 / r a / rs rs cs/Wi cs(Te/mi)1/2
• ITER 1 / r 1000
• Simplest model GK ions, adiabatic electrons
• Ion Temperature Gradient (ITG) driven
  turbulence
• Relevant scales krrs, kqrs lt 1
• 3D field solver min 4 pts/wavelength
• Nr (2/p) 1/r 640, Nq 4 1/r 4000, Nj
  Nq/q 3000
• Phase space resolution 300 markers / cell
  
• or Nm8, Nv//32
• Needs to be assessed by numerical convergence
  studies
• Time until statistical steady-state 2000 a /
  cs
• Needs to be assessed by series of (long) runs
  with independent initial conditions

Nm11
Fourier mode in filter
71
ITER global GK turbulence simulation
requirements, with field-aligned Fourier
20 x 106
11

• Field solver 600 x 4000 x 3000 7.2 x109 grid
  points.
• CPU cost scaling (1/r)3
• 2000 x 109 particles (Lagrange-PIC)
• CPU cost scaling (1/r)3
• 2000 x 109 grid points (Eulerian / semi-Lagrange
  field-aligned grid)
• CPU cost scaling (1/r)3
• Number of time steps, assuming Wi Dt 1, Nsteps
  t/Dt 2000 (a/cs) Wi 2000 (1/r)
  2x106
• CPU cost scaling (1/r)
• The computational cost scales as (1/r)4 and is
  prohibitive for ITER-size parameters
• Improvements must be found !

(1/r)2
6
(1/r)2
6
(1/r)2
Dt cs/a 0.1
2 x 104
(1/r)2
(1/r)0
Feasible on top-end existing HPCs ()
have been
()needs scalable massively parallel algorithms
72
A question of size, benchmarks, and

• GK turbulence simulation codes are often run at
  the limit of available computing power available.
  Verification vs analytical theories exist only in
  simplified cases. Hence, cross-code comparisons
  (benchmarks) are essential.

• CYCLONE base case A.M. Dimits et al, PoP 7
  (2000) 969
• Local (flux tube) results (in approximate
  geometry) compared with global code results. Good
  agreement.
• Crucial question how finite size simulations
  (i.e. global) tend to infinite system size (i.e.
  flux tube) cases?

Cited 306 times as of Aug.31, 2009 ISI Web of
Science
73
A question of size, benchmarks, and
A controversy had risen regarding the size
scaling of turbulent-driven transport
J. Candy, R.E. Waltz, W. Dorland, PoP 11 (2004)
L25
GTC global Lagrange-PIC GYRO global Euler
Dimits flux-tube infinite system size limit
74
A question of size, benchmarks, and
1/r a/rs ITG
Inconsistent geometrical approximation
smaller physical system size

• Numerically converged results (Np, Dt) within
  symbol size
• Agreement between global and flux tube results
  for CYCLONE DIII-D parameters Dimits PoP 2000
  is purely coincidental!

X. Lapillonne et al., PoP 16, 032308 (2009)
75
A question of size, benchmarks, and
consistent exact geometry
consistent approximate geometry
inconsistent approximate geometry
Fig.6 X. Lapillonne et al., PoP 16, 032308
(2009)

• Most of the difference can be attributed to
  inconsistency of the geometrical approximation

76
A question of size, benchmarks, and

• Stiffness of anomalous transport implies very
  large sensitivity in the vicinity of marginality!

77
ORB5 a Lagrangian, PIF code Particle-In-Finite
Element

• Based on GK conservative formulation of TS Hahm
• Axisymmetric tokamak geometry, link to MHD
  equilibrium
• Fully global, f evolution is unconstrained
• allows for profile relaxation, avalanches,
  bursts, etc
• Control variates only dff-f0 is sampled with
  markers
• adaptive / enhanced control variates
• Importance sampling optimized loading
• Finite element 3D cubic B-spline for dn, df
  fields
• Straight-field-line coordinates with
  Field-Aligned Fourier Filtering
• Order gain in computing cost scaling with system
  size (r)
• Noise-control with sources, zonal flow preserving
• Long, statistically converged simulations above
  marginal stability

78
Zonal Flows self regulate turbulence
Color contours of perturbed electrostatic
potential
hot core
S. Jolliet, et al., CPC 2007
large scale Zonal Flows gt breakup of turbulent
eddies gt self-regulation of turbulence
cold edge
turbulent eddies gt radial transport

• Direct Numerical Simulation of Ion Temperature
  Gradient driven turbulence
• Resolution 128x448x320, N 83 million, WiDt40
• 6 h on 1024 processors IBM BG/L_at_EPFL
• No sources final state has (undamped) Zonal
  Flows and decaying turbulence level

79
Interaction of zonal flows with turbulence can
been observed in many other physical systems
Zonal flow velocity profile of Jupiters high
atmosphere
Zonal flow velocity profile in a cylindrical
plasma
NASA Cassini 2000
Idomura, et al., JCP 2007
Quasi-2D electrostatic plasma turbulence
Quasi-2D atmospheric turbulence
80
Coordinate singularity near axis
Pseudo-cartesian coordinates (x,h) grid
Polar-like coordinates (s,q) grid

• Equation of motion for q is singular on axis

• Equations of motion for x and h are regular

81
Outline

• 1. Turbulence and transport, GyroKinetic theory
• 2. Numerical schemes to solve GK equations
• 3. Critical issue no.1 numerical noise
• 4. Critical issue no.2 geometry, anisotropy
• 5. Critical issue no.3 long, (quasi-)steady
  state simulations
• Chaotic nature of transport processes
• Filamentation of phase space in dissipationless
  systems
• Necessity of finite dissipation for a
  steady-state
• Active noise control algorithms
• Entropy production, entropy balance
• 6. Collisions in PIC
• 7. High Performance Computing issues
• 8. Outlook and conclusions

82
4. Critical issue no.3 long, (quasi-)
steady-state simulations

• Turbulent processes have a chaotic nature
• Predictions made by direct numerical simulations
  (DNS) of turbulence have a statistical aspect
• Aim obtain long enough, statistically converged,
  steady-state simulations
• Some dissipation must be present in the system in
  order to get a statistical steady-state Krommes
  Hu, PoP 1 (1994) 3211, Krommes PoP 6 (1999)
  1477
• Dissipationless (collisionless) systems lead to
  indefinitely thin filamentation of f in phase
  space. In PIC, this leads to an indefinite
  increase in the fluctuation entropy and decay of
  the Signal/Noise ratio (though lower order
  moments of f, e.g. the energy flux, may seem to
  have converged entropy paradox)

83
Turbulent processes have a chaotic nature

• Convergence studies with increasing Np
• Nonlinear ITG, sourceless, collisionless, white
  noise initialisation

S. Jolliet, et al, CPC 177 (2007) 409
84
Turbulent processes have a chaotic nature

• Convergence studies with increasing Np
• Nonlinear ITG, sourceless, collisionless, mode
  initialisation

Signature of chaos finite time after which two
neighboring initial conditions diverge from each
other
S. Jolliet, et al, CPC 177 (2007) 409
85
Turbulent processes have a chaotic nature

• Convergence studies with increasing Np
• Nonlinear ITG, with sources (see further)

B.F. McMillan, et al, PoP 15 (2008) 062306
How to define a good (converged???) estimate of
the flux? Time-averages?
86
Turbulent processes have a chaotic nature

• Varying initialisations
• Nonlinear ITG, with sources (see further)

B.F. McMillan, et al, PoP 15 (2008) 062306
Temporal moving average over 500 a/cs of
normalised heat diffusivity
87
Collisionless, dissipationless filamentation

• Core fusion plasmas are very weakly collisional
  (lmfp,// gtgt device size)
• Issue for computations achieving numerical
  convergence is impossible after a finite time
  (resolution should increase indefinitely) in
  PIC, this manifests itself as a decaying
  Signal/Noise ratio
• Finite dissipation is necessary

v
Spatial scales bounded by Larmor radius and
device size gt Direct Numerical Simulations (DNS)
are possible, but Collisionless (Landau)
wave-particle interactions lead to filamentation
in phase space This filamentation proceeds
indefinitely
x
88
Statistical noise in PIC simulations
N markers in phase space
Nm Fourier modes in the field representation
Coeff of order 1 dependent on FLR and spline
functions

• In collisionless, dissipationless nonlinear
  simulations, the noise increases indefinitely
  because ltw2gt increases indefinitely
• Noise-control put a limit (bound) on the
  increase of w2

89
Long, statistically converged simulations

• Add sources/sinks in order to
• Maintain average gradients above marginal
• Avoid accumulation of error (noise)
• Preserving the physics, in particular Zonal Flows

Krook avoid noise accumulation
correction to preserve undamped Zonal Flows
heating
overbar orbit average
lt.gt surface average
90
Zonal Flow phase space structure

• Contours of df in ZF damping test using the
  global Eulerian code GT5D Idomura et al., CPC
  2008

Coherent structure due to neoclassical
polarisation effect
Filament structures produced by ballistic modes
91
Zonal-flow preserving corrected Krook
ORB5
Damped GAM
ltvExBgt a.u.
McMillan, et al., PoP 2008
t Wci-1

• The undamped residual Rosenbluth-Hinton is
  preserved

92
A long nonlinear simulation with sources
Heat transport
Signal / noise

• Numerical noise is controlled over very long
  times, using the noise-control Krook operator

93
With/without Krook as noise control
McMillan, et al., PoP 2008

• When S/N is lost, there is a decrease in computed
  transport level

94
Noise destroys the field structure
Using a heating operator noise-dominated
Using a Krook operator noise-controlled
McMillan, et al., PoP 2008

• Contours of non-zonal perturbed potential f - ltfgt

95
Checking quasi-state entropy balance

• The fluctuation entropy, i.e. difference between
  macroscopic and microscopic entropies, is defined
  as Watanabe, Sugama, PoP 9 (2002) 3659

• From the gyrokinetic equation, we obtain the
  following entropy balance relation

entropy production from profile gradients
dissipative noise control term (lt0)
Dfield particles to field perturbation (usually
small)
Dheat contribution of heating operator (usually
small)
96
Checking quasi-stateness entropy balance

• ORB5 simulation, with heating, no noise-control

97
Checking quasi-stateness entropy balance

• ORB5 simulation, with heating, with noise-control

98
Bursts, avalanches role of Zonal Flows

• A substantial fraction of heat transport is
  carried by bursts

GAMs
In/out propagating bursts avalanches
dltfgt/ds
Quasi-steady Zones
radial coordinate

• NB similar avalanches are seen on the heat flux,
  ZF shearing rate, and gradT

99
Outline

• 1. Turbulence and transport, GyroKinetic theory
• 2. Numerical schemes to solve GK equations
• 3. Critical issue no.1 numerical noise
• 4. Critical issue no.2 geometry, anisotropy
• 5. Critical issue no.3 long, (quasi-)steady
  state simulations
• 6. Collisions in PIC
• Langevin (random walk) approach
• Weight spreading, noise
• 2-weights scheme
• 7. High Performance Computing issues
• 8. Outlook and conclusions

100
6. Collisions in PIC

• Even though core fusion plasmas are weakly
  collisional, finite collisionality has effects on
  e.g.
• Reducing Trapped Electron Mode (TEM) drive
• Damping the residual Zonal Flows (ZF)
• In the absence of turbulence, in curved magnetic
  confinement systems, they lead to the so-called
  neo-classical transport (NC)
• NC transport include heat, particle, and momentum
  transport (e.g. the bootstrap current effect,
  very important for advanced tokamak scenarios)
• In transport barriers (regions where
  turbulence-induced transport is suppressed) NC
  transport remains
• For turbulence simulations, they bring physics
  first-principle based dissipation mechanisms
  (which are, as we have just seen, essential in
  order to obtain steady-state simulations)

101
Collisions Langevin approach

• The effect of particle collisions is a diffusion
  and a drag
• For example, electron-ion collisions can be
  modeled with the Lorentz operator, which is a
  limiting case of the more general Landau operator

Langevin random walk kick
102
Collisions Langevin approach

• In df schemes, the particle weights, wp, have to
  be interpreted statistically
• It is like an additional dimension in phase space
• 2-weight scheme G. Hu and J.A. Krommes, PoP 1
  (1994) 863
• The Langevin approach leads to weight
  spreading, hence to increased particle
  statistical noise
• Ways have been devised to limit the detrimental
  effects of weight spreading, e.g. an evolving
  background scheme
• S. Brunner, E. Valeo, J.A. Krommes, PoP 6 (1999)
  4504

103
Collisions Langevin weight spreading
Relax the individual weights w toward their
averaged values W
104
Collisions neoclassical effects

• Considering collisions alone (no turbulence!)

R/a2.8, a/LT0, a/Ln5
105
Collisions Trapped Electron Mode

• electron-ion collisions decrease TEM instability
  drive

106
Collisions (ion-ion) zonal flow damping

• Ion-ion collisions damp the residual zonal flows

107
Collisions (i-i) and turbulent transport
Z. Lin, et al, Phys. Rev. Lett. 83 (1999) 3645,
PoP 7 (2000) 1857

• Increased collisionality -gt more damping of Zonal
  Flows -gt increased turbulence level -gt increased
  transport level

108
Outline

• 1. Turbulence and transport, GyroKinetic theory
• gyrokinetic modelling of fusion plasmas
• 2. Numerical schemes to solve GK equations
• Lagrange (PIC), Euler, Semi-Lagrange
• 3. Critical issue no.1 numerical noise
• Statistical interpretation of PIC
• Control Variates
• Loading issues
• 4. Critical issue no.2 geometry, anisotropy
• Straight-field-line coordinates
• Field-aligned Fourier filter
• Benefits for noise problem
• Benefits for field solver problem
• Beware of simplified geometrical models
• 5. Critical issue no.3 long, (quasi-)steady
  state simulations
• Chaotic nature of transport processes
• Problem of dissipationless systems filamentation
  of phase space
• Necessity of finite dissipation for a
  steady-state
• Active noise control algorithms

109
Outline

• 1. Turbulence and transport, GyroKinetic theory
• 2. Numerical schemes to solve GK equations
• 3. Critical issue no.1 numerical noise
• 4. Critical issue no.2 geometry, anisotropy
• 5. Critical issue no.3 long, (quasi-)steady
  state simulations
• 6. Collisions in PIC
• 7. High Performance Computing issues
• Massive parallelism
• Parallel scalability
• 8. Outlook and conclusions

110
7. High Performance Computing

• Even with the recent improvements in algorithms,
  direct numerical simulation of turbulence using
  gyrokinetic models is very demanding in computer
  resources
• The evolution of top-end HPC to massively
  parallel platforms requires the corresponding
  development of parallelizable and scalable
  algorithms
• www.top500.org the top of the list consists of
  platforms having several 100,000s cores
• In the not-so-distant future millions of cores!
• Not so easy to tap the immense cpu power of such
  architectures
• Requires re-thinking of algorithms (re-factoring
  of codes)

111
Parallelization scheme ORB5
Aim solve very large problems gt huge nr of
procs gt huge nr of domains !?!!!?
Domain decomposition (toroidal direction)
MPI move
particles
Domain Cloning field quantities (density,
potential) are replicated
Nc 1 100
MPI global sum
fields
Npes NdNc
Nd 100 1000
112
Massive parallelism

• Lagrange ORB5 code. Strong scaling with 8x108
  particles (left). Weak scaling with 8x105
  particles/processor (right).
• IBM BG/L in co-processor mode
• 6.4 G particles enough to simulate ITG
  turbulence with adiabatic electrons in the whole
  core of ITER
• Thanks to a DEISA project (Distributed European
  Infrastructure for Supercomputing Applications),
  scalability was demonstrated up to 32768
  processors.

113
Massive parallelism - Scalability

• BG/P (Idris) and XT5 (Monte Rosa, CSCS)

114

• Thanks T. Stitt, CSCS DEISA-DECI

115
8. Conclusion and outlook

• Algorithm improvements have resulted in orders of
  magnitude increase in performance
• These have been based on a good understanding of
  both physical and computational aspects
• The capability to control error (noise)
  accumulation has been shown, allowing for long,
  quasi-steady state simulations
• Scalability with system size is close to ideal
• Parallel scalability demonstrated up to 32768
  processors and promising
• Global, gyrokinetic Direct Numerical Simulation
  of turbulence in the ITER plasma core appears
  feasible in a near future
• at least for Ion-Temperature-Gradient driven
  turbulence
• Work is underway to include and study more
  physical effects
• ElectroMagnetic, Collisions, Fast ions (e.g.
  fusion alphas),
• Code development and cross comparisons will
  continue
• HPC issues (massive //, ) will continue to
  dominate our agenda
• Require human resource development on a long term
  perspective

116
Fusion
Credit NASA, ESA and J. M. Apellániz (IAA, Spain)

• in space (a nursery of stars NGC 6357)

117
Controlled fusion on earth magnetic confinement
Source EFDA-JET
JET has produced 16MW of fusion power
Source EFDA and J.B. Lister CRPP

• Figure of merit density x Temperature x energy
  confinement time (nTtE )
• Magnetic confinement tE1s, n1020m-3, T10keV.
  ITER nTtE3.4 atm.s

118
Semi-Lagrange approach

• Main motivation get the better of two worlds
• avoid the statistical noise problem inherent of
  Lagrange-PIC-MC methods
• avoid the CFL condition inherent of Euler methods

• GYSELA code
• 5D gyrokinetic
• Operator splitting
• Leapfrog scheme, with averaging of integer and
  half-integer timesteps

fconst along orbit
Orbit, backward integrated from grid point
Interpolation of f _at_ foot of orbit
119
Lagrange (ORB5) vs ½-Lagrange (GYSELA)
Alternative approaches bring useful cross-checks
of numerical and physical results
Heat diffusivity
Ion temperature gradient
120
Outline

• 1. Introduction turbulence and transport
• gyrokinetic modelling of fusion plasmas
• 2. Numerical schemes to solve GK equations
• Lagrange, Euler, Semi-Lagrange
• 3. Critical issues
• anisotropy, role of zonal flows, numerical error
  (noise)
• 4. Quasi-steady state simulations
• sources/sinks, cascades, avalanches
• 5. High Performance Computing
• massive parallelism, scalability
• 6. Outlook and conclusions