Gyrokinetic Simulations Using Particles - PowerPoint PPT Presentation

1 / 113
About This Presentation
Title:

Gyrokinetic Simulations Using Particles

Description:

(1) Centre de Recherches en Physique des Plasmas, Association Euratom - Suisse, ... Properties of GK-Maxwell equations (in the collisionless, sourceless case) ... – PowerPoint PPT presentation

Number of Views:64
Avg rating:3.0/5.0
Slides: 114
Provided by: CRPP3
Category:

less

Transcript and Presenter's Notes

Title: 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
Write a Comment
User Comments (0)
About PowerShow.com