Kinetic Simulation of Sheared Flows in Tokamaks

1
Kinetic Simulation of Sheared Flows in Tokamaks

• J.A. Heikkinen2, S.J. Janhunen1, T.P. Kiviniemi1,
  S.Leerink1, M. Nora1, and F. Ogando1,3,
• 1 Teknillinen Korkeakoulu (TKK), Finland
• 2 Valtion Teknillinen Tutkimuskeskus (VTT),
  Finland
• 3 Universidad Nacional de Educación a Distancia
  (UNED), Spain

2
Outline

• Background
• Gyrokinetic full f approach for toroidal fusion
  plasmas.
• Examples of Vlasov codes.
• Description of a PIC full f code.
• Testing
• Comparison to neoclassical theory.
• Linear and nonlinear benchmarks of unstable
  modes.
• Influence of noise on results.
• Transport simulations in toroidal plasmas.
• Future challenges of full f

3
Section IBackground
4
Turbulence and zonal flows

• Self-organization and zonal flow generation with
  turbulence is important not only in fusion
  plasmas but also in other fluids like planetary
  atmoshpheres and solar plasma
• Transport transients and related transport
  barriers are important, e.g., for ITER fusion
  reactor design
• Empirical evidence from a number of tokamak
  facilities indicates a complex scaling law for
  transport transients
• The scaling laws for transport barriers in
  tokamaks are today to a large extent not
  understood by physics principles

5
Sheared flows and turbulence appear in toroidal
plasma simulation
6
Kinetic simulation required

• All trials up to date around the world to see the
  spontaneous confinement transition in an
  experimental-like first principles plasma
  simulation have failed so far.
• Reason to this may be that not all relevant
  physics like kinetic details of particle motion
  have been included in the simulation codes, lack
  of self-consistency with background evolution, or
  model simplifications like using fluid codes have
  been too crude.

7
Section IIGyrokinetic full f approach for
simulation of toroidal plasmas
8
Full f simulation of tokamak plasmas

• In full f simulation, one solves for the whole
  particle distribution function in phase space and
  in time either non-perturbatively, with
  gyrokinetics or drift-kinetically
• Pioneering non-perturbative 5D particle-in-cell
  (PIC) simulation of magnetically confined
  toroidal plasmas (unrealistic me/mi)
• Cheng Okuda, 78 Sydora, 86 LeBrun, 93,
  Kishimoto, 94
• Eulerian and semi-Lagrangian Vlasov 3D, 4D, and
  5D simulation
• Cheng Knorr, 76 Manfredi, 96
  Sonnendruecker, 99 Xu, 06 Scott, 05
  Grandgirard, 06
• Gyrokinetic 5D particle simulation
• Furnish, 99 Heikkinen, 00, 03 Chang, 03.

9
Calculation flux in particle codes
Calculation of forces from fields and velocity
Acceleration and increment of velocity
Computation of electric field. Magnetic is given.
Initial step with ?0
Displacements and new positions. Boundary
conditions.
Calculation of density. Current profile fixed.
Resolution of Poisson equation for the
electrostatic potential.
10
Delta f vs. full f

• Delta f calculates perturbations from an assumed
  background distribution f0.
• Powerful for small f-f0
• Linear mode analysis
• Snapshot transport analysis
• Path-breaking global transport studies for large
  toroidal installations

• Full f calculates the whole particle
  distribution.
• Fitting processes that perturbate strongly the
  particle distribution
• Strong transient or long time scale transport in
  core or edge plasmas
• Strong particle/energy sources
• Full neoclassical equilibrium
• Edge plasma (sheaths, wall losses, recycling,
  separatrix, flows)
• Large MHD (sawteeth, ELMs)

11
Delta f vs. full f

• Automatic extension of Delta f to full f
  calculation may be possible
• Non-conservative effects collisions, sources
• Loading optimization
• Automatic increase of the number of markers when
  required
• Few particles (10-100) per cell are needed for
  good results with small f-f0.

• Full f can be realized both in PIC and Vlasov
  (Eulerian, semi-Lagrangian)
• Vlasov simulation is free of noise and has
  controlled numerical resolution in phase space.
• PIC is numerically flexible and straightforward
  to implement
• Needs many particles per cell (1000) or fine
  grid discretization of velocity space for an
  acceptable noise level or accuracy.

12
Section IIIExamples of full f Vlasov codes
13
Full f 5D gyrokinetic code GYSELA Grandgirard,
06 (CEA, LPMIA-Univ, IRMA-Univ, LSIIT-Illkirch)

• Semi-Lagrangian GYSELA code
• - fixed grid, follows trajectories backwards,
• - global code
• - 5D GAM and ITG Cyclone tests successful

14
Full f 5D gyrokinetic code TEMPEST Xu 06 (LLNL,
Calif-Univ, GA, LBNL, PPPL)

• - Based on modified gyrokinetics valid for large
  long-wavelength and small short-wavelength
  variations (Qin, 06)
• - Fixed grid, equations solved via a
  Method-of-Lines approach and an implicit
  backward-differencing scheme using iteration
  advances the system in time
• - Developed for circular core or divertor edge
  geometry
• - 4D neoclassical tests successful. 5D drift wave
  and ITG benchmarking ongoing

15
Section IVThe ELMFIRE code an example of a full
f PIC code
16
The ELMFIRE group
Founded in 2000 at Euratom-Tekes Contributors
from Finland Spain Holland Main
affiliations VTT TKK ... but also ... CSC Åbo
Akademi UNED (Spain)
17
ELMFIRE code

• Full f nonlinear gyrokinetic particle-in-cell
  approach for global plasma simulation (present
  version electrostatic).
• Magnetic coordinates (?,?,?) Boozer 81.
• Guiding-center Hamiltonian White Chance, 84.
• Gyrokinetics is based on Krylov-Boholiubov
  averaging method in description of FLR effects
  (P. Sosenko, 01).
• Adiabatic or kinetic electrons with impurities.
• Parallelized using MPI with very good
  scalability.
• Based on free software PETSc and GSL for math
  calc.

18
ELMFIRE full f features

• An initial canonical distribution function avoids
  the onset of unphysical large scale ExB flows
  (Heikkinen, 01)
• Direct implicit ion polarization (DIP) and
  electron acceleration (DEP) sampling of
  coefficients in the gyrokinetic equation
• Quasi-ballooning coordinates to solve the
  gyrokinetic Poisson equation

• Versatile heat (RF, NBI, Ohmic) sources and
  particle sources/ recycling
• Full binary collision operator

19
Poisson equation

• W.W. Lee proposed standard model with
  polarization drift included in equation operator.
• Ion density evaluated from ion motion without
  polarization drift
• P. Sosenko proposes including polarization in the
  ion density.
• Ion density evaluated from ion motion with
  polarization drift. Circular gyro-orbits.

20
Implementation into ELMFIRE

• Solve ? by isolating ion polarization drift
  contribution to density.
• That contribution is calculated implicitely every
  timestep using also previous values of ?.
• The gyroaveraged electric field is interpolated
  from grid potential values for the ion
  polarization drift.

Larmor circle
ith subparticle cloud
y
k
i
?py
?px0
xp,yp
?px-
?px
?py0
ds
k
?py-
x
pth point on the Larmor orbit of the kth ion
gyro-center of the kth ion
21
Section VBenchmarking of ELMFIRE
22
Testing ELMFIRE

• Comparison to neoclassical theory in the presence
  of turbulence.
• Frequency of GAM and Rosenbluth residual.
• Neoclassical radial electric field.
• Comparison to other codes has been done in the
  Cyclone Base cases.
• Linear growth of unstable modes and their phase.
• Nonlinear saturation of transport in both
  adiabatic and kinetic-electron case.
• Comparison to experimental results.
• Collaboration with IOFFE Institute and St.
  Petersburg Polytechnic working with the FT-2
  tokamak.

23
Available resources

• Gyrokinetic full f computation is very demanding
  computationally. Parallel computation is a need.
• CSC (The Finnish IT Center for Science) provides
  shared use of high-end parallel computers.

• IBM eServer cluster 1600. 512 processors with
  2.2TFlops, 384GB RAM and High Performance Switch
  communication.
• Cluster of 768 AMD OpteronTM processors up to
  3.2TFlops, 1600GB RAM, Infiniband network.
• Cray XT4 (Hood) 70TFlop, 70TB RAM HP ProLiant
  Supercluster, 10.6 Tflop, 100 TB

24
Geodesic Acoustic Modes

• Neoclassical theory predicts GAM frequency and
  Rosenbluth residual.
• Results show good wide agreement with theory.
• Simulations done on a plasma annulus.
• R0.3-0.9 m, a0.08 m, B0.6-2.45 T, q1.28-2.91,
  Ti90-360 eV, ni5.11019m-3 (r/a0.75)

25
Neoclassical radial electric field

• Neoclassical radial electric field is well
  followed in conventional (L-mode) turbulent
  simulations both in radius and in time

• R1.1 m, a0.08 m, B2.1 T, I22 kA, parabolic
  ion heated n,Ti,e profiles (r0.04m).

26
Neoclassical benchmarking

• Parameters and boundary conditions
• FT-2 like parameters, R00.55 m and a0.08 m,
• Itot 55 kA, T,n,j (1-(r/a)2)? ,
• n051019 1/m3
• T0300 eV in high T case, T0120 eV in low T
• no heating, no loop voltage
• relaxing profiles, cooling by recycling on
  outer radius

27
Neoclassical benchmarking
Number of particles must be sufficient as Er may
depend on noise Er radial dependence fairly well
predicted by standard neoclassical theory.
However, Reynolds stress and poloidal Mach number
can be important.
28
Linear growth of unstable modes

• Test based on adiabatic Cyclone Base Case (Dimits
  PoP '00)
• Red points from ELMFIRE, blue line fit from
  article.
• Figures show growth rates and typical time
  evolution for a mode with k??i0.3

29
Linear growth of unstable modes

• Test based on adiabatic Cyclone Base Case (Dimits
  PoP '00)
• Red points from ELMFIRE, blue line fit from
  article.
• Figures show growth rates and typical time
  evolution for a mode with k ?? i0.3

Region of linear growth
30
Linear growth of unstable modes

• Test based on kinetic Cyclone Base Case (Chen NF
  '03)
• Filled circles and squares from ELMFIRE
• Dashed line fit for the growth rate ? from Chen
  NF 03.

31
Evolution of thermal conductivity

• Evolution of ?i is studied with nonlinear runs of
  Cyclone Base.
• Measured at ra/2 (q1.4). Using kinetic
  electrons. R/LT10. Weak collisionality
  T(a/2)2000 eV, n51017 m-3.
• Convergence requires a large number of particles
  per cell.

32
Section VIInfluence of noise on results
33
Influence of noise on results

• Particle simulation not only covers physical
  density fluctuations, but it produces undesirable
  noise with fluxes that perturbate the solution.
• Associated diffusivity can be estimated from the
  radial particle shift during decorrelation time.
• Physical radial ion heat conductivity can be
  estimated from mixing-length estimate of the
  physical level of fluctuations, being also
  proportional to T3/2.

34
Effects on calculated conductivity

• Image shows influence of strong collisionality
  and potential averaging on ion radial heat
  conductivity.
• Collisionless cases show residual noise
  conductivity.
• Noise is filtered out by averaging potential over
  flux surface.
• So far noise is reduced by brute force ...
  higher N!

Scaled Cyclone Base Case with kinetic electrons
T100 eV, n4.51019 m-3
35
Contribution of noise to the heat flux
The convective noise flux prediction gives a
fairly good estimate of the unphysical ion heat
conductivity in the simulations
36
From noise to turbulence spectrum
Wave number spectrum from different time
instants demonstrates how one moves from white
noise to physical spectrum in modes.
37
Problematic levels of noise
Fluctuations _at_ ra/2

• Figure show exceptionally bad case regarding
  noise effects.
• It is a kinetic cyclone base case with scaled
  parameters and low T100eV and n4.51017 m-3.
• Density fluctuations remain almost constant in
  time.
• Regression shows almost perfect N-1/2 scaling,
  indicating that results are dominated by noise.

38
but not always problematic...

• In FT-2, fluctuation levels are much higher
  (10-40) than in scaled Cyclone Base Case (1).
• So high perturbation level warrants the use of a
  full f scheme.
• Image shows density fluctuations relative to flux
  surface average.
• Relative importance of noise values can be seen
  in videos of both cases.

39
Probability distribution functions

• Fluctuation levels can be represented by the PDF
  graphs.
• PDFs show fluctuation distributions over a middle
  flux surface.
• Density values are averaged over time.
• Turbulence in FT-2 takes fluctuations up to 40
  in start (up) and 15 after 50µs (down).

40
Section VIITransport simulations
41
Case 1 under study LH heated FT-2

• Parameters from 100 kW LH heated 22 kA FT-2
  tokamak.
• Case shows a rapid growth of potential and
  electric field and reduction of thermal
  conductivity interpreted as ITB formation.
• Top figure diagram (r,t) of flux surface average
  of potential.
• Bottom figure ion diffusivity at mid radius.

42
Evolution of profiles

• Strong ExB flow shear is created at r0.05 m,
  where a knee-point in Ti profile is found

43
Calculated spectra

• Parameters of the FT-2 case, devised to cause the
  formation of Internal Transport Barrier.
• B2.2T, I22kA, n3.51019m-3, Ti250eV, Te300eV

S(k) at r5.1cm
? vs. time
S(k) at r7.5cm
44
Case 2 under study LH heated FT-2

• Heating phase for 100 kW LH heated 22 kA FT-2
  tokamak (O8 impurities included).

45
Evolution of large scale fluctuations

• Density fluctuations plots show the formation and
  further destruction of macroscopic structures

46
Evolution of diffusivity

• Both particle diffusivity and heat conductivity
  drop drastically when poloidal flow shear
  destroys the turbulent structures
• The figures show values from the middle radius

47
Case 3 under study Ohmic FT-2

• Reproduce the FT-2 reflectometer signal I(t)
  with ELMFIRE
• I(t)? w(r,?) dn(r,?,t) r dr d?
• W(r, ? ) Weighting function calculated by beam
  tracing code using exp. data
• dn(r,?,t) Density fluctuations simulated by
  ELMFIRE code
• Results
• Poloidal velocity calculated by the shift of the
  power spectrum of I(t) is in reasonable agreement
  with the poloidal velocity measured by the
  Doppler reflectometer for an Ohmic discharge
• Width of the power spectrum is too narrow
  compared to the experimental results.

48
Frequency broadening is too narrow in the
simulation
49
Section VIIIFuture challenges
50
Unresolved issues in full f

• How to ensure a sufficient number of particles in
  all grid cells when density varies strongly in
  global PIC simulation
• Adaptation of good grid resolution for strongly
  varying f in Vlasov simulations
• Full collision operator in Vlasov simulations
• SOL plasma simulation with sheath boundaries
• Gyrokinetics with strong fluctuations

51
Scalability
IB-new Sepeli with InfiniBand and optimized
parallelization IB Sepeli with InfiniBand LAM
Sepeli with Gigabit Ethernet IBM-new IBM with
optimized parellelization
52
Resource scenario
53
Conclusions

• Successful 5D gyrokinetic full f PIC and Vlasov
  simulations of tokamak electrostatic turbulence
  and transport for core plasma.
• Careful benchmarking of the codes is performed in
  appropriate limits for the turbulence saturation
  and neoclassical characteristics.
• Linear and nonlinear benchmarking.
• Vlasov approach has less noise and better control
  of resolution PIC code is more flexible to
  implement
• Most urgently needed for edge plasma simulations
  further development of nonlinear terms in
  gyrokinetics may be needed

ACKNOWLEDGEMENTS
This project receives funding from the European
Commission
CSC The Finnish IT Center for Science for its
computing facilities