INTRODUCTION TO GYROKINETIC AND FLUIDSIMULATIONS OF PLASMA TURBULENCEAND OPPORTUNITES FOR ADVANCED FUSIONSIMULATIONS
G.W. Hammett, Princeton Plasma Physics Labw3.pppl.gov/∼hammettFusion Simulation Project Workshop
San Diego, Sept. 17, 2002

Thanks to Bill Nevins &the Plasma Microturbulence Projectfor many vugraphs. In particular see:
http://www.isofs.info/nevins.pdf
http://fusion.gat.com/theory/pmp
and others working on Braginskii
2-fluid simulations of edge turbulence.
(Collabs. at Univ. Alberta, UCLA,UCI, Univ. Colorado, Dartmouth, Garching,
General Atomics, LLNL, Univ. Maryland,MIT, Princeton PPPL, Univ. Texas)A complete description of a plasma
is given by the particle distribution function F s(~x, ~v , t), the density of particles
at (near) position ~x with velocity ~v and time t, for species s (with charge q s andmass ms).
The charge density and current needed for Maxwell’s equations to determine
the electric and magnetic fields is then:σ(~x, t) =XsZ~j(~x, t) = X qs d3v~v Fs(~x, ~v , t)

3qs d vFs(~x, ~v , t)sZ

Fs is determined by the Vlasov-Boltzmann equation~  ∂F∂F~v × B∂Fqs ·~+ ~v ·+= Collisions + sources + sinks ≈ 0
E +∂t∂~x msc∂v
where sources + sinks includes radiation cooling of electrons, ionization andrecombination changes of ion charge state, etc.

Equivalent particle approach
Discrete particle density representation (combined with smoothing and“particle-in-cell” techniques):
Fs =d~xi= ~vidtX

i=1,Nwi (t)δ(~x − ~xi(t)δ(~v − ~vi(t)~ xi, t) d~vi~vi × B(~qi ~xi, t) +=E(~dtmicplus Monte Carlo treatment of collisions, sources and sinks.
Both “continuum” F and particle descriptions are equivalent (in the limit ofa large number of particles, typical fusion particle density ∼ 10 14/cm3) andare“Exact”, but both include an excessive range of time and space scales.
Most plasma phenomena of interest are slow compared to the electron and ion
gyrofrequencies (∼ 10 11 Hz and ∼ 108 Hz).

The Nonlinear Gyrokinetic EquationGuiding center distribution function F s(~x, ~v , t) = F0s (ψ, W ) + F0s (ψ, W )qsφ̃/Ts +h̃s(~x, W, µ, t) = equilibrium + fluctuating components, where the energy W =mvk2 + µB, the first adiabatic invariant µ = mv ⊥2 /B, and∂ h̃s ˜
∂F0s ∂ χ̃+ ~v χ + vk b̂ + ~vd · ∇h̃s = −~v˜χ · ∇F0s − qs∂t
∂W ∂t+ Collisions + Sources + Sink

where b̂ points in the direction of the equilibrium magnetic field, ~v d is the curvature and grad B drift, Ω s is the gyrofrequency, and the ExB drift is combined withtransport along perturbed magnetic fields lines and the perturbed ∇B drift as:vk ! J1(γ) mv⊥2 B̃kχ̃ = J0(γ) φ̃ − Ãk +cγe Bc~v˜χ = b̂ × ∇χ̃B

J0 & J1 are Bessel functions with γ = k ⊥ v⊥/Ωs, and the fields are from
∂F0s0 ≈ 4π qs d v qsφ̃+ J0(γ)h̃ss∂W
4π X Z 32∇ Ãk = −qs d vvk J0(γ)h̃sc sJ1(γ)
B̃k4π X Z 3=− 2h̃sd vmv⊥2B
B sγ
XZ3



~˜ = −∇φ̃ − 1 ∂ Ãk b̂
Ec ∂t~ =B~ 0 + ∇Ãk × b̂ + B̃kb̂
B
In a full-torus simulation where plasma variations must be kept
1 ZJ0(k⊥ v⊥/Ωs)φ → hφi(~x) =d~ρφ(~x + ρ~)
2π.

Candy/Waltz movies available at:
http://web.gat.com/comp/parallel/gyro gallery.html
and other movies can be found from various links starting at:
http://fusion.gat.com/theory/pmp

ThePlasma
Microturbulence
Project supports a 2x2 matrix
of codes (geometry x algorithm),
each type of code is tuned tooptimize in various regimes andso are optimized to study certain
types of problems.
Codes using flux-tube geometry(shown here) take advantage ofshort decorrelation lengths of theturbulence perpendicular to magnetic field lines. Multiple copies ofa flux-tube pasted together represent a toroidal annulus.

We Support a 2x2 Matrix ofPlasma Turbulence Simulation Codes•Continuum

PIC

Flux TubeGS2SUMMIT

GlobalGYRO

GTC

Why both Continuum and Particle-in-Cell (PIC)?
– Cross-check on algorithms– Continuum currently most developed (already has kinetic e’s , B )– PIC may ultimately be more efficient•

If we can do Global simulations, why bother with Flux Tubes?
– Electron-scale ( e, e=c/ pe) physics (ETG modes, etc.)– Turbulence on multiple space scales (ITG+TEM, TEM+ETG, ITG+TEM+ETG, …)– Efficient parameter scans

June 5, 2002

Plasma Microturbulence Project3

Current ‘state-of-the-art’(similar performance achieved in Continuum codes)•Spatial Resolution
Plasma turbulence is quasi-2-D–

Resolution requirement along B–fielddetermined by equilibrium structure– Resolution across B–field determinedby microstructure of the turbulence.
⇒ ~ 64×(a/ρi)2 ~ 2×108 grid points tosimulate ion-scale turbulence atburning-plasma scale in a global code
– Require ~ 8 particles / spatial grid point⇒ ~ 1.6×109 particles for global ionturbulence simulation at ignition scale– ~ 600 bytes/particle⇒ 1 terabyte of RAM•

Temporal Resolution
Studies of turbulent fluctuations–

Characteristic turbulence time-scale⇒ cs/a ~ 1 µs (10 time steps)– Correlation time >> oscillation period⇒ τc ~ 100× cs/a ~ 100 µs(103 time steps)– Many τc’s required⇒ Tsimulation ~ few ms(5×104 time steps)– 4×10-9 sec/particle-timestep(this has been achieved)⇒ ~90 hours of IBM-SP time/run⇒ This resolution is achievable* Heroic (but within our time allocation)(Such simulations have been performed, see T.S. Hahm, Z. Lin, APS/DPP 2001)•Simulations including electrons and B (short space & time scales) are not
yet practical at the burning-plasma scale with a global code

5/23/02

W.M. Nevins, "Turbulence & Transport"7

Major Computational andApplied Mathematical Challenges•Continuum kernels solve an advection/diffusion equation on a 5-D grid
– Linear algebra and sparse matrix solves (LAPAC, UMFPAC, BLAS)– Distributed array redistribution algorithms (we have developed or own)•

Particle-in-Cell kernels advance particles in a 5-D phase space– Efficient “gather/scatter” algorithms which avoid cache conflicts and provide
random access to field quantities on 3-D grid••Continuum and Particle-in-Cell kernels perform elliptic solves on 3-D grids(often mixing Fourier techniques with direct numerical solves)
Other Issues:––––

Portability between computational platforms
Characterizing and improving computational efficiency
Distributed code development
Expanding our user base

5/23/02

W.M. Nevins, "Turbulence & Transport"

19

Continuum / Eulerian Codes
Flux-tube / thinannulus

Full-torus or thin annulus

Particle-in-Cell/Lagrangian Codes
Flux-tube

Full-torus

All now use field-line following coordinate systems, ∆x⊥/∆x|| ~ ρi/L ~ 10-1-10-3GS2 (Dorland, U. Md., Gyro (Candy-Waltz GA)
Kotschenreuther)

GTC (Z. Lin et.al.
Summit(LLNL, U. Co, PPPL, UCI)
UCLA)

⊥ Pseudo-spectral
linear & nonlinear.

Delta-f algorithm reduces particle
noise.

Toroidal pseudo-spectral

5th-6th order upwind τ grid
Recent hybrid electron algorithm:|| 2cd order finite-diff. to avoid 1/v||
fluid with kinetic electron closure.
(slight upwind)collisions w/ direct sparse
solver (UMFPACK)
Linear: fully implicit(elegant algorithm)
Nonlinear: 2cd order
Adams-BashforthHigh accuracy explicit 4thorder Runge-Kutta

Elliptic solvers easyin Fourier space

Elliptic solvers with non-uniform coefficients solved by
combination of Fourier, iterative, and direct matrix solution

Leap-frog / Predictor-corrector

Fast time scales hiding in E & B fields: is there a partially-implicit iterativealgorithm that can help?

Recommendations (I)Strengthening PMP Support to Integrated Modeling(1) Improve the fidelity and performance of Plasma
Microturbulence Project codes(2) Validate these codes against experiment(3) Expand the user base of the PMP codes(4) Initiate the development of a kinetic edge
turbulence simulation code.

5/23/02

W.M. Nevins, "Turbulence & Transport"

13

CORE TURBULENT TRANSPORT STILL IMPORTANT• Provides most of temperature gradient: 20 keV center → 1-4 keV near-edge.
Effects of shaping, density peakedness, rotation, impurities, T i /Te ?
• Detailed experimental comparisons possible, fluctuation diagnostics.
• Are internal transport barriers possible at reactor scales? P threshold? Torque?
Controllable?
• Electron-scale transport controls advanced reactor performance?

BUT EDGE TURBULENCE CRITICAL• H-mode pedestal (edge transport barrier) greatest source of uncertainty forreactor predictions.
• Will divertor melt/erode? Need ELM simulation.
• Edge very complicated: Separatrix & divertor geometry matters. Bootstrap
current important, second stability regime. Half of power radiated, intense
neutral recycling.
• High and low collisionality regimes. Present edge codes are collisional fluids, need kinetic extensions.

3-D Fluid Simulations ofPlasma Edge Turbulence
BOUT (X.Q. Xu, )•Braginskii — collisional, two fluid
electromagnetic equations•

Realistic ×-point geometry(open and closed flux surfaces)•

BOUT is being applied to DIII-D,
C-Mod, NSTX, …

•

There is LOTS of edge fluctuation data!

⇒ An Excellent opportunity forcode validation
5/23/02

W.M. Nevins, "Turbulence & Transport"

15

More info:
Plasma Microturbulence Project (PMP):
http://fusion.gat.com/theory/pmp
Nevins presentation on PMP to ISOFS May 2002:
http://www.isofs.info/nevins.pdf
GS2 (Dorland Univ. Md.):
http://gk.umd.edu/GS2/info.html
Useful 2-page gyrokinetic summary:
http://gk.umd.edu/GS2/gs2 back.ps
GTC (Lin PPPL UCI):
http://w3.pppl.gov/∼zlin/visualization/
Gyro (Candy/Waltz GA):
http://web.gat.com/comp/parallel/gyro.html
Summit (LLNL/UCLA/U. Co.):
http://www.nersc.gov/scidac/summit