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 v H ≈ v A = B/ √ 4 πn i m i . • 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 • Shear Alfv´ en waves in non-uniform equilibria exhibit a continu- Z ous spectrum: “local” plasma oscillations with frequency contin- r uously changing throughout the plasma ϑ • In toroidal geometry, the poloidal-symmetry breaking due to the R toroidal field variation on a given magnetic flux surface cause R a 0 different poloidal harmonics to be coupled: frequency “gap” ap- pears in the Alfv´ en continuum ϕ • Discrete, global MHD modes (Toroidal Alfv´ en Eigenmodes, or 1 ω TAE’s) can exist in the gaps of the shear-Alfv´ en frequency spec- m=2 0.8 trum. TAE’s are marginally stable MHD modes and can be easily driven unstable by the resonance with energetic particles 0.6 m=1 • Use reduced MHD equations expanded up to O ( ǫ 3 ), with ǫ ≡ 0.4 a/R 0 , a and R 0 the minor radius and the major radius of the 0.2 m=1 torus, respectively, to keep toroidal effects in the model m=2 0 0 0.2 0.4 0.6 0.8 1 r 4

• Hybrid MHD-kinetic models – Energetic particle density is typically much smaller than the bulk plasma density – Ordering: n H T H ≈ O ( ǫ 3 ) , ≈ O ( ǫ − 2 ) n i T i – Thus, the following ordering for the ratio of the energetic to bulk ion β ( β ≡ 8 πP 0 /B 2 0 is the ratio between the plasma kinetic and the magnetic pressures) follows: β H ≈ O ( ǫ ) β i – 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: Discretized form of F ( t, Z ): ∂t + dZ i ∂F ∂F � dZ ′ F ( t, Z ′ ) δ ( Z − Z ′ ) ≈ F ( t, Z ) ≡ ∂Z i = 0 , dt N l =1 ∆ l F ( t, Z l ) δ ( Z − Z l ) . � Equations of motion: dZ i dt = . . . . Phase-space grid points Z l ( 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). dZ i   d ∆ l  ∂ Volume elements ∆ l ( t ) evolve according to: dt = ∆ l ( t ) .   ∂Z i dt  t,Z l ( t ) Often it could be convenient to evolve only the perturbed part δF of the distribution function: = ⇒ F ( t, Z ) = F 0 ( t, Z l ) + δF ( t, Z l ) 6

4 Particle-in-cell versus Finite-size-particle Plasma condition n 0 λ 3 D ≫ 1 (collective interaction dominate over collisions) implies a huge number of simulation particles. Even assuming n s λ 3 D ≈ 10, typically ( L eq equilibrium length): D ( L eq /λ D ) 3 ≈ 10 13 . N part ≈ n s L 3 eq = n s λ 3 Violation of plasma condition n 0 λ 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 – L c – 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 n 0 L 3 = c ≫ 1 (with L c ≫ λ D ) 7

• Finite-size-particle (FSP) Finite-size-particles (charge clouds): n s ( x ) = l ∆ l δ ( x − x l ) − → l ∆ l S ( x − x l ) � � δ (x) S(x) L s x x The spatial characteristic width of the cloud λ D ≪ L s ≪ L eq restricts the maximum spatial resolution attainable in the simulation (assuming L eq /L s ≈ 100, n s λ 3 D ≈ 10): s ( L eq /L s ) 3 ≈ 10 7 . n s λ 3 → n s L 3 N part ≈ n s L 3 D ≫ 1 − s ≫ 1 , L s plays the role of L c 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: O PIC ≈ f ( N harm ) + n FT × N harm × N cell + n int × N part , M PIC ≈ m harm × N harm + m cell × N cell + m part × N part , N harm : number of Fourier harmonics retained in the simulation; f ( N harm ): operations for the solution of the field solver; n FT : number of operations needed to compute each addendum in the Fourier transform; N cell : number of cells of the spatial grid; n int : operations for the field interpolation; N part : number of simulation particles; m harm , m cell and m part : 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: O FSP ≈ f ( N harm ) + n FT × N harm × N part , M FSP ≈ m harm × N harm + m part × N part . Typically, f ( N harm ) negligible in comparison with terms ∝ N part ; for PIC codes, N ppc ≡ N part /N cell ≈ n 0 L 3 c ≫ 1: as far as n FT × N harm ≫ n int 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, or weak nonlinear coupling (nonlinear mode spectrum restricted to a limited number of signif- icant harmonics): in such a few-harmonic framework, the condition n FT × N harm ≫ n int can be violated or, at least, significantly weakened; 2. schemes of parallelization (distributed-memory architectures): domain decomposition ( d.d. ) versus particle decomposition ( p.d. ) i procs =4 i procs =3 i procs =2 i procs =1 i procs =1 i procs =2 i procs =3 i procs =4 10

Different portions of the physical domain are assigned to Domain decomposition, PIC different processors, together with the particles that reside on them. d.d. ≈ f ( N harm )+ 1 O PIC ( n FT × N harm × N cell + n int × N part ) , n proc d.d. ≈ m harm × N harm + 1 M PIC ( m cell × N cell + m part × N part ) . n proc Advantages: almost linear scaling of the attainable physical-space resolution (more precisely, the maximum number of spatial cells) with the number of processors. i procs =1 i procs =2 i procs =3 i procs =4 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 of a PIC code very complicate. 11