The Hybrid MHD-Gyrokinetic Code HMGC G. Vlad * * Associazione - - PowerPoint PPT Presentation

The Hybrid MHD-Gyrokinetic Code HMGC

• G. Vlad*

*Associazione Euratom-ENEA sulla Fusione,

C.R. Frascati, Rome, Italy

in collaboration with S. Briguglio*, G. Fogaccia*, F. Zonca* and B. Di Martino†

†Second University of Naples, Naples, Italy

III Convegno Nazionale su “La Fisica del Plasma in Italia” L’Aquila, 20-22 Maggio 2002

Electronic version: http://fusfis.frascati.enea.it/˜vlad/Miscellanea/slides L’Aquila 2002.pdf

1

1 Outline

• Introduction
• The model
• Numerics

– Fluid section – Gyrokinetic section: particle simulations, Particle-in-cell (PIC) vs. Finite-size-particle (FSP) – Parallelization: Domain vs. Particle decomposition – Parallel architectures: Distributed Memory, Shared Memory, Hierarchical Distributed- Shared Memory

• Examples

2

2 Introduction

• The Hybrid MHD-Gyrokinetic Code (HMGC) has been developped at the C.R. Frascati,

ENEA laboratory in the frame of thermonuclear fusion research

• Recent experimental devices are approaching the so called ignition condition: fusion α-particles

are confined in the toroidal (Tokamak) plasma and sustain the burning plasma

• Confinement properties of the energetic (α) particles are crucial in obtaining good performances

in reactor relevant regimes

• Fusion α-particles are born mainly in the plasma centre, and the corresponding radial profile

is peaked: their pressure gradient is a free-energy source that can destabilize waves which resonantly interact with the periodic motion of the energetic particles

• Energetic

(hot) ions (α-particles) in plasmas close to ignition conditions have vH ≈ vA = B/√4πnimi.

• Interaction between energetic particles and shear Alfv´

en waves is likely to occur

• Shear Alfv´

en waves = ⇒ Magnetohydrodynamic (MHD) model

• Energetic particles (wave-particle interaction) =

⇒ Kinetic model

• Confinement properties of energetic particles =

⇒ Nonlinear model

3

3 The model

• Bulk plasma: described by Magnetohydrodynamic (MHD) equations

R a ϑ r R Z ϕ

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 ω r m=2 m=1 m=1 m=2

• Shear Alfv´

en waves in non-uniform equilibria exhibit a continu-

• us spectrum: “local” plasma oscillations with frequency contin-

uously changing throughout the plasma

• In toroidal geometry, the poloidal-symmetry breaking due to the

toroidal field variation on a given magnetic flux surface cause different poloidal harmonics to be coupled: frequency “gap” ap- pears in the Alfv´ en continuum

• Discrete, global MHD modes (Toroidal Alfv´

en Eigenmodes, or TAE’s) can exist in the gaps of the shear-Alfv´ en frequency spec-

• trum. TAE’s are marginally stable MHD modes and can be easily

driven unstable by the resonance with energetic particles

• Use reduced MHD equations expanded up to O(ǫ3), with ǫ ≡

a/R0, a and R0 the minor radius and the major radius of the torus, respectively, to keep toroidal effects in the model

4

• Hybrid MHD-kinetic models

– Energetic particle density is typically much smaller than the bulk plasma density – Ordering: nH ni ≈ O(ǫ3) , TH Ti ≈ O(ǫ−2) – Thus, the following ordering for the ratio of the energetic to bulk ion β (β ≡ 8πP0/B2

0 is

the ratio between the plasma kinetic and the magnetic pressures) follows: βH βi ≈ O(ǫ) – It can be shown that, making use of the above ordering, the MHD momentum equation is modified by a term which represent the perpendicular component of the divergence of the energetic-particle stress tensor ΠH – Energetic-particle stress tensor obtained by solving Vlasov equation

5

• Particle simulation: gyrokinetic model

Direct solution of the equation describing the time evolution of the particle distribution function F(t, Z) for collisionless plasmas: Vlasov equation: ∂F ∂t + dZi dt ∂F ∂Zi = 0 , Equations of motion: dZi dt = . . . . Discretized form of F(t, Z): F(t, Z) ≡

• dZ′F(t, Z′)δ(Z − Z′) ≈

N

• l=1 ∆lF(t, Zl)δ(Z − Zl) .

Phase-space grid points Zl(t) evolve according to eqs. of motion: numerical particles. Gyrocenter coordinates Z ≡ (R, µ, v, θ): R is the gyrocenter position, v parallel (to B) velocity, µ magnetic moment (exactly conserved in this coordinate system), θ the gyrophase (does not appear explicitly). Volume elements ∆l(t) evolve according to: d∆l dt = ∆l(t)

   ∂

∂Zi dZi dt

  

t,Zl(t)

. Often it could be convenient to evolve only the perturbed part δF of the distribution function: = ⇒ F(t, Z) = F0(t, Zl) + δF(t, Zl)

6

4 Particle-in-cell versus Finite-size-particle Plasma condition n0λ3

D ≫ 1 (collective interaction dominate over collisions) implies a huge number

• f simulation particles. Even assuming nsλ3

D ≈ 10, typically (Leq equilibrium length):

Npart ≈ nsL3

eq = nsλ3 D (Leq/λD)3 ≈ 1013 .

Violation of plasma condition n0λ3

D ≪ 1: system too collisional, short range interactions between

particles dominate over the long range ones. Particle-in-cell (PIC)

• 1. Electromagnetic fields computed at the points of a discrete spatial

grid

• 2. Interpolation of the e.m. fields at the (continuous) particle posi-

tions to compute the forces and perform particle pushing

• 3. Pressure contribution of energetic particle calculated at the grid

points to close the equations = ⇒ Short-range interactions are then cut off for mutual distances shorter than the typical spacing – Lc – between grid points, whilst the relevant long-range interactions are not significantly affected. = ⇒ PIC particle ensemble behaves as a plasma under the much more relaxed condition n0L3

c ≫ 1

(with Lc ≫ λD)

7

• Finite-size-particle (FSP)

Finite-size-particles (charge clouds): ns(x) =

• l ∆lδ(x − xl)

− →

• l ∆lS(x − xl)

δ(x) x Ls S(x) x

The spatial characteristic width of the cloud λD ≪ Ls ≪ Leq restricts the maximum spatial resolution attainable in the simulation (assuming Leq/Ls ≈ 100, nsλ3

D ≈ 10):

nsλ3

D ≫ 1 −

→ nsL3

s ≫ 1 ,

Npart ≈ nsL3

s (Leq/Ls)3 ≈ 107 .

Ls plays the role of Lc

8

5 Computational loads and Parallelization Assume that field solver uses Fourier transform to solve the MHD equations. Serial code, number of operations (O) per time step and memory (M) required: PIC: OPIC ≈ f (Nharm) + nFT × Nharm × Ncell + nint × Npart , M PIC ≈ mharm × Nharm + mcell × Ncell + mpart × Npart ,

Nharm: number of Fourier harmonics retained in the simulation; f (Nharm): operations for the solution of the field solver; nFT: number of operations needed to compute each addendum in the Fourier transform; Ncell: number of cells of the spatial grid; nint: operations for the field interpolation; Npart: number of simulation particles; mharm, mcell and mpart: memory needed to store, respectively, a single harmonic of the complete set of Fourier-space fields, the real-space fields at each grid point and the phase-space coordinates for each particle.

FSP: OFSP ≈ f (Nharm) + nFT × Nharm × Npart , M FSP ≈ mharm × Nharm + mpart × Npart. Typically, f (Nharm) negligible in comparison with terms ∝ Npart; for PIC codes, Nppc ≡ Npart/Ncell ≈ n0L3

c ≫ 1: as far as nFT × Nharm ≫ nint the gridless FSP method is more

expensive than the PIC one, without presenting any significant advantage in terms of memory requests.

9

Two distinct reasons could however justify a different trend:

• 1. Interest in simplified simulations in which only very few modes are evolved: linear simulations,
• r weak nonlinear coupling (nonlinear mode spectrum restricted to a limited number of signif-

icant harmonics): in such a few-harmonic framework, the condition nFT × Nharm ≫ nint can be violated or, at least, significantly weakened;

• 2. schemes of parallelization (distributed-memory architectures):

domain decomposition (d.d.) versus particle decomposition (p.d.)

iprocs=1 iprocs=3 iprocs=2 iprocs=4

iprocs=1 iprocs=3 iprocs=2 iprocs=4

10

Domain decomposition, PIC

iprocs=1 iprocs=3 iprocs=2 iprocs=4

Different portions of the physical domain are assigned to different processors, together with the particles that reside

• n them.

OPIC

d.d. ≈ f (Nharm)+ 1

nproc (nFT × Nharm × Ncell + nint × Npart) , M PIC

d.d. ≈ mharm×Nharm+ 1

nproc (mcell × Ncell + mpart × Npart) . Advantages: almost linear scaling of the attainable physical-space resolution (more precisely, the maximum number of spatial cells) with the number of processors. Disadvantages: particle migration from one portion of the grid to another, possible severe load-balancing problems = ⇒ dynamical redistribution of grid and particle quanti- ties is required, which makes the parallel implementation

• f a PIC code very complicate.

11

Particle decomposition, PIC

iprocs=1 iprocs=3 iprocs=2 iprocs=4

Statically distributing the particle population among pro- cessors, while replicating the data relative to grid quanti-

• ties. Before updating the electromagnetic fields, at each

time step, partial contributions to particle pressure coming from different portions of the population must be summed together (reduction). OPIC

p.d. ≈ f (Nharm) + nFT × Nharm × Ncell + nint × Npart

nproc , M PIC

p.d. ≈ mharm × Nharm + mcell × Ncell + mpart × Npart

nproc . Advantages: load balancing is automatically enforced; par- allelization is, in principle, almost straightforward. It is very efficient if computational load related to particles dominates, for each processor, the one related to the grid (nproc< ∼ Nppc). Disadvantages: grid calculations do not take advantage, with regard both to the number of operations and the mem-

• ry requests, from such a parallelization: each processor has

to handle the whole spatial domain. Even neglecting effi- ciency problems, high spatial-resolution levels are limited by the single node RAM.

12

Particle decomposition, FPS

iprocs=1 iprocs=3 iprocs=2 iprocs=4

bottle-necks in efficiency and performance associated to grid quantities induces one to by-pass the introduction of a spa- tial grid, so resorting to the gridless FSP simulation: OFSP

p.d. ≈ f (Nharm) + nFT × Nharm × Npart

nproc , M FSP

p.d.

≈ mharm × Nharm + mpart × Npart nproc . Advantages: Fourier transforms are distributed among the

• processors. High spatial resolution can be obtained. Mas-

sively parallel simulations can yield significant benefits as far as the number of modes, Nharm, retained in the simula- tion is relatively small, in spite of the high mode numbers considered (high spatial resolution). Disadvantages: few-harmonic limitation.

13

6 HMGC Parallel architectures

• The HMGC code exists in a PIC version and in a gridless FSP version.
• Parallel implementations include:

– Distributed Memory (IBM SP, cluster of workstations), – Shared Memory (Symmetric Multiprocessor Systems, SMPs), – Hierarchical distributed-shared memory multiprocessor architectures (IBM SP, cluster of SMPs).

Level INTER-NODE Language HPF Strategy Particle Decomposition Level INTRA-NODE Language OpenMP Phase Particle pushing Pressure updating Variant Version Strategy Particle decomposition Particle decomposition Domain decomposition Critical Auxiliary array p aux Sorting Selective sorting – v1 v2a v2b 14

7 Results Fluid nonlinearities: saturation of a TAE (ω = ω0 ≈ 0.33 ωA)

10-14 10-12 10-10 10-8 10-6 10-4 100 200 300 400 500 WTOT m,n ω A t (1,0) (1,1) (2,1) (3,2)

0.30 0.32 0.34 0.36 0.38 0.40 0.66 0.68 0.70 0.72 0.74 ω/ω A s linear phase non-linear phase

Volume integrated energy (magnetic plus kinetic) for different Fourier components (m, n) vs. time for a non-linear simulation of an unstable driven TAE. The q-profile has a parabolic radial depen- dence with q(0) = 1.1 and q(a) = 1.9. The inverse aspect ratio is ǫ = 0.075, the density is constant ˆ ̺ = ˆ ̺0 and the resistivity corresponds to S−1 = 10−5. Blow-up of the Alfv´ en continuum. The continuous spectra obtained in the linear limit and at the beginning of the non-linear phase are compared.

15

Kinetic nonlinearities: gap-mode saturation (ω = ω0 ≈ 0.33 ωA)

6 4 2 200 400 600 800 ωA t ln A (t) 0.8 0.7 0.6 0.5 2π 4π r/a 0.8 0.7 0.6 0.5 r/a Ψ2,1 2π 4π Ψ2,1 a) b)

Time evolution of the mode amplitude A(t) for a perturbative non-linear simulation with γD = 0.01 ωA, βH(0) = 0.08. Orbit in the plane (Ψ2,1, r) (Ψm,n ≡ ωrt − mϑ + nϕ), for a test particle, in the time intervals 0 < ωAt < 284 (a) (linear growth) and 264 < ωAt < 560 (b) (non-linear saturation). The Ψ2,1 axis is mapped onto the interval 0 ≤ Ψ2,1 < 4π. The particle is initially passing, but becomes trapped as the mode reaches a certain amplitude.

16

Transition from Kinetic TAE to EPM (Energetic Particle Mode):

0.00 0.04 0.08 0.12 0.0 0.2 0.4 0.6 0.8 0.005 0.015 0.025 0.035 βH γ /ω A ω r/ω A ω r KTAE ω r EPM γ KTAE γ EPM a) 0.0 0.2 0.4 0.6 0.8 1.0 r/a b)

β-threshold vs. toroidal mode number n:

0 ,1 0 ,2 0 ,3 0 ,0 1 0 ,0 2 0 ,0 3 0 ,0 4 0 ,0 5

β

H

γτ

A n=8 n=4 n=1 0.00 0.01 0.02 0.03 5 1 0 1 5 2 0

n βH th

P I C F S P

17

Nonlinear EPM saturation generating shear flows (n = 8, monotonic q(r) profile, q(0) = 1.1, q(a) = 1.9, Npart ≈ 16.7 × 106): see movie: http://fusfis.frascati.enea.it/˜vlad/Miscellanea/EPM MOVIES/n8 9 imirr1 13 zonal 3x4.mov

18

Deeply hollow q profile

• Deeply hollow q profile (q(0) ≈ 5, qmin = 2.1, q(a) ≈ 5, profile (a)), βH(0) = 2.5%.
• ωgap/ωA,r=0 = 1/(2q(r)
• ρ/ρr=0): assume first a radially constant thermal-plasma density ρ

⇒ radially constant Alfv´ en velocity (such an assumption will be removed later).

see movie: http://fusfis.frascati.enea.it/˜vlad/Miscellanea/IAEA-Goteborg/n4 JET 7.mov

• After a transient initial phase, a mode localized

around the maximum β′

H emerges at r ≈ 0.35a,

with frequency well inside the continuum. We can identify this mode as an Energetic Particle contin- uum Mode (EPM).

• Its saturation takes places because of a strong

(convective) outward radial displacement of the energetic ions.

• As such a displacement takes place, the lo-

cal drive is reduced due to the flattening of the energetic-ion density profile. The drive is no longer able to overcome the strong continuum damping at the original frequency.

• The maximum of the power spectrum migrates

towards the gap (in order to minimize the contin- uum damping), but it also moves outwards, fol- lowing the displaced source, in order to maximize the drive.

• The mode reaches the gap and it localizes

around the zero-shear, qmin surface (r ≈ 0.5a).

19

R a ϑ r R Z ϕ

Low-β tokamak ordering (β ≡ 8πP0/B2

0 is the

ratio between the plasma kinetic and the mag- netic pressures): v⊥ vA ≈ B⊥ Bϕ ≈ B/B · ∇ ∇⊥ ≈ O(ǫ) , vϕ vA ≈ ∇ · v⊥ vA/a ≈ ∇(RBϕ) Bϕ ≈ O(ǫ2) , ∂ ∂t ≈ vA R . Here, a cylindrical-coordinate system (R, Z, ϕ) has been used, and the subscript ⊥ denotes com- ponents perpendicular to ∇ϕ. The magnetic field can be written as B =

• F0 + ˜

F

• ∇ϕ + R0∇ψ × ∇ϕ + O(ǫ3Bϕ)

where ψ is the poloidal magnetic flux function, F0 = R0B0, B0 is the vacuum (toroidal) magnetic field at R = R0, and ˜ F ≈ O(ǫ2F0) is given, at the leading order, by equilibrium corrections.

20

• Reduced MHD equations:

∂ψ ∂t = − cR2 R0B0 ∇ψ × ∇ϕ · ∇φ − c R0 ∂φ ∂ϕ + η c2 4π∆∗ψ + O(ǫ4vABϕ) , ˆ ̺

  D

Dt − 2c R0B0 ∂φ ∂Z

  ∇2

⊥φ + ∇ˆ

̺ ·

  D

Dt − c R0B0 ∂φ ∂Z

  ∇φ =

− B0 4πcB · ∇∆∗ψ − B0 cR0 ∇ ·

• R2 (∇P + ∇ · ΠH) × ∇ϕ
• + O(ǫ4̺v2

ABϕ

a2c ) , DP Dt = O(ǫ4vAB2

ϕ

a ) , with η the plasma resistivity, v⊥ = − cR2 R0B0 ∇φ × ∇ϕ + O(ǫ3vA) , ˆ ̺ = R2 R2 ̺ , D Dt = ∂ ∂t + v⊥ · ∇ , ∇2

⊥ ≡ 1

R ∂ ∂RR ∂ ∂R + ∂2 ∂Z2 , ∆∗ψ = R2∇ ·

 ∇ψ

R2

  = R ∂

∂R

  1

R ∂ψ ∂R

  + ∂2ψ

∂Z2 .

21

• Particle simulation target: obtaining from the numerical plasma the same behaviour of the

physical one It is impossible to simulate, with today numerical resources, the number of particles of real plasmas Consider the Debye length λD and the Larmor radius ρL: λD =

 

T 4πe2n

 

2

, ρL = c(mT)2 eB Consider mutually interacting macroparticles: ns = nf M , M ≫ 1 , es = Mef , ms = Mmf , vs = vf = ⇒ Ts ∝ msv2

s = Mmfv2 f ∝ MTf

= ⇒ λD,s = λD,f , ρL,s = ρL,f Plasma condition n0λ3

D ≫ 1 (collective interaction dominate over collisions) implies a huge number

• f simulation particles. Even assuming nsλ3

D ≈ 10, typically (Leq equilibrium length):

Npart ≈ nsL3

eq = nsλ3 D (Leq/λD)3 ≈ 1013 .

Violation of plasma condition n0λ3

D ≪ 1: system too collisional, short range interactions between

particles dominate over the long range ones.

22

The pressure tensor is: ΠH (t, x) = 1 m2

H

• d6ZDzc→ZF H (t, R, M, U) ×

  ΩHM

mH I + ˆ bˆ b

  U2 − ΩHM

mH

      δ (x − R) ,

Equation of motion for energetic particles: dR dt = Uˆ b + eH mHΩH ˆ b × ∇φ − U mHΩH ˆ b × ∇a +

   M

mH + U ΩH

 U + a

mH

     ˆ

b × ∇ ln B, dM dt = 0, dU dt = 1 mH ˆ b·

        eH

ΩH

 U + a

mH

  ∇φ + M

mH ∇a

  

× ∇ ln B + eH mHΩH ∇a×∇φ

   − ΩHM

mH ˆ b · ∇ ln B.

23