• No results found

The Hermite-Fourier spectral method for solving the

N/A
N/A
Protected

Academic year: 2021

Share "The Hermite-Fourier spectral method for solving the"

Copied!
48
0
0

Loading.... (view fulltext now)

Full text

(1)

IN , SECOND DEGREE PROJECT ELECTROPHYSICS 120 CREDITS

CYCLE

STOCKHOLM SWEDEN 2016 ,

The Hermite-Fourier spectral method for solving the

Vlasov-Maxwell system of equations

JURIS VENCELS

KTH ROYAL INSTITUTE OF TECHNOLOGY

(2)

Sammanfattning

Denna avhandling fokuserar p˚ a f¨ orb¨ attring av Hermite-Fourier spektral metod f¨ or att l¨ osa kinetiska plasma problem. I den f¨ orsta delen av avhandlingen en ny dynamiskt adaptiva tekniker f¨ or att ¨ andra antalet av Hermite l¨ agen och Hermite grund presenteras. F¨ orkonditionering av problemet och anv¨ andning av vetenskaplig toolkit PETSc diskuteras. Tekniken f¨ or att ¨ andra antalet l¨ agen och prekonditionering tros minska ber¨ akningstid, medan en f¨ or¨ andring av grunden f¨ orb¨ attrar robusthet och numerisk konvergens av metoden. Den andra delen av avhandlingen fokuserar p˚ a parallellisering strategier och prestandaanalys av koden implementeras i programmeringsspr˚ aket Fortran.

Abstract

This thesis focuses on the improvement of the Hermite-Fourier spectral method

for solving kinetic plasma problems. In the first part of the thesis a novel

dynamically adaptive techniques for changing the number of Hermite modes

and Hermite basis are presented. Preconditioning of the problem and use of

high-end scientific toolkit PETSc are discussed. The technique of changing the

number of modes and preconditioning are believed to reduce computational

time, while a change of basis improves robustness and numerical convergence of

the method. The second part of the thesis focuses on paralellization strategies

and performance analysis of the code implemented in the Fortran programming

language.

(3)

Contents

1 Introduction to Computer Simulations of Plasma 7

1.1 Numerical Methods for Plasma Physics . . . . 8

1.2 Scope of Thesis . . . . 9

2 Literature Survey and Theoretical Background 10 2.1 The Hermite-Fourier Spectral Method . . . . 10

2.2 Velocity Expansion in Hermite Basis . . . . 10

2.3 Spatial Expansion in Fourier Basis . . . . 12

2.4 Filamentation Problem and Collisional Operator . . . . 14

2.5 Implicit Time Discretization . . . . 15

2.6 Transformation to Real Space . . . . 15

2.7 Extension to the Multi-Dimensional Case . . . . 16

3 Development of the Hermite-Fourier Method 19 3.1 Dynamically Adaptive Number of Fluid Moments . . . . 19

3.1.1 Dynamically Adaptive Change of the Number of Hermite Modes . . . . 20

3.1.2 Results . . . . 20

3.2 Dynamically Adaptive Change of Hermite Basis . . . . 22

3.2.1 Basis Transformation Method . . . . 23

3.2.2 Results . . . . 27

3.3 Non-linear Solvers and Preconditioners . . . . 31

3.3.1 Description of PETSc . . . . 32

3.3.2 Construction of Preconditioner . . . . 32

3.3.3 Results . . . . 34

4 Paralellization of the Method and Performance Optimization 37 4.1 Profiling of the 1D1V and 2D3V Implementations . . . . 37

4.1.1 Allinea MAP . . . . 37

4.1.2 Profiling Results . . . . 38

4.2 Paralellization with MPI on Distributed-Memory Systems . . . . 40

4.3 Paralellization with OpenMP on Shared-Memory Systems . . . . 41

4.4 Scaling Test Results . . . . 41

5 Discussion and Conclusions 43 List of Figures 1 First six Hermite basis functions Ψ n (ξ). . . . 11

2 Adaptive number of Hermite modes for Landau damping and two stream instability tests. . . . . 21

3 Savings in terms of time and Hermite modes for technique of adaptive number of Hermite modes . . . . 22

4 Basis shift convergence example. . . . 25

5 Basis scaling convergence example. . . . 25

6 Distribution function error for different scaling and shift param-

eters. . . . 26

(4)

7 Two stream instability test for simulations with and without dy-

namically adaptive change of basis. . . . 28

8 Space-phase plot for two-stream instability with and without dy- namically adaptive change of basis. . . . 29

9 Non-linear Landau damping tests with and without dynamically adaptive change of basis. . . . 30

10 Space-phase plot for non-linear Landau damping with and with- out dynamically adaptive change of basis. . . . 31

11 Dispersion relation of the ion acoustic wave. . . . 35

12 Electric and magnetic fields for Whistler instability. . . . . 36

13 Profiling results for 1D1V implementation . . . . 39

14 Profiling results for 2D3V implementation . . . . 40

15 Scaling test for 1D1V implementation . . . . 42

16 Scaling test for 2D3V implementation . . . . 42

List of Tables 1 Initial parameters for the Landau damping with ions and the two stream instability simulations. . . . 20

2 Simulation parameters for adaptive change of basis. . . . 27

3 Simulation parameters for ion acoustic wave with preconditioned PETSc solver . . . . 34

4 Comparison of simulation time for preconditioned and not pre- conditioned 1D1V ion acoustic wave simulation for different time steps. . . . . 35

5 Simulation parameters for whistler instability with preconditioned PETSc solver . . . . 36

6 Comparison of simulation time for preconditioned and not pre- conditioned 2D3V whistler instability simulation for different time steps. . . . . 36

7 Data access latencies for different kinds of memories. . . . 38

8 Two stream instability simulation parameters for performance profiling. . . . 39

9 Whistler instability simulation parameters for performance pro-

filing. . . . . 40

(5)

Nomenclature

Acronyms

API Application Program Interface

CN Crank-Nicolson

DRAM Dynamic random-access memory

GMRES Generalized Minimal Residual

MKL Hermite-Fourier

HPC High Performance Computing

JFNK Jacobian-Free Newton-Krylov

MHD Magnetohydrodynamics

MKL Math Kernel Library

OpenMP Open Multi-Processing PDE Partial Differential Equation

PETSc Portable, Extensible Toolkit for Scientific Computation

PIC Particle-in-Cell

Notation

α Velocity scaling parameter

C m (x, t) Expansion coefficient in Hermite basis

C mn (t) Expansion coefficient in Hermite and Fourier basis

e Elementary charge

f (x, v, t) Distribution function

H m m-th order Hermite polynimial

L Length of spatial domain

λ D Debye length

m Particle’s mass

N The highest Fourier mode

n Plasma density

N F = 2N + 1 The total number of Fourier modes N H The total number of Hermite modes

q Particle’s charge

(6)

s Species

u Velocity shift parameter

v t Thermal velocity

ω Plasma frequency

ξ Transformed velocity (scaled and shifted) Normalization 1D

E = ˆ m E e

e

v

te

ω

pe

Electric field f = ˆ f v n

te

ref

Distribution function, where n ref is a reference plasma density

ˆ

m s = m s /m e Mass normalized by electron mass ˆ

q s = q s /e Charge normalized by elementary charge

ˆ t = t ω pe Time normalized by inverse electron plasma frequency ˆ

v = v/v te Velocity normalized by the electron thermal velocity ˆ

x = x/λ D Coordinates normalized by electron Debye length Normalization 2D3V

B = ˆ B B

0

Magnetic field normalized by reference magnetic field E = ˆ c B E

0

Electric field

ˆ t = t ω pe Time normalized by inverse electron plasma frequency ˆ

v = v/c Velocity normalized by the speed of light ˆ

x = x ω

pe

c Coordinates normalized by electron inertial length

(7)

1 Introduction to Computer Simulations of Plasma

The name “plasma” in its physical context was first used in 1927 by Ameri- can chemist and physicist Langmuir to describe an ionized gas [1]. Since then plasma physics has become a fundamental research discipline that is extremely important for everyday technologies and industrial applications.

Fusion devices for energy generation, propulsion systems for space explo- ration, techniques for plasma etching and sputtering to fabricate integrated cir- cuits, and security in space environment are the main plasma-related research topics today.

Although laboratory, space, and astrophysical plasmas often consist of sim- ple components like electrons, protons and heavy nuclei, the mass difference encompassing at least three orders of magnitude drives different plasma dy- namics occurring over a large range of time and spatial scales. Because of this multi-scale nature, numerical modelling of plasmas is extremely challenging.

Depending on the state, plasma can be described with different models. For example, the Klimontovich equation [2]

∂N s (x, v, t)

∂t + v · ∂N s (x, v, t)

∂x + e s m s

[E(x, t) + v × B(x, t)] · ∂N s (x, v, t)

∂v = 0, (1) where N s (x, v, t) = P N

0

n=1 δ[x − x n (t)]δ[v − v n (t)] is microscopic distribution function, N 0 the total number of particles for species s. This equation describes behaviour of individual particles and is exact description of plasma. In practice it is rarely used in computer simulations because of high computation costs and problem sizes.

Alternatively, plasma can be described by macroscopic velocity moments which lead to a so-called “fluid” model. For example, Magnetohydrodynam- ics (MHD) [3] describe plasma with three macroscopic quantities - density, bulk velocity and pressure. This description is advantageous from computation view- point since it reduces the problem size and complexity, but only plasmas with Maxwellian velocity distribution (strong collisionally) can be studied. The MHD model also has a number of constraints in length and time scales.

Use of higher velocity moments can reduce these time and spatial limitations of MHD, but derivation of generalized higher-moment fluid equations is a very challenging task and it is not always clear how many moments are needed for a given problem.

The Vlasov equation, the most popular equation in plasma physics, is con- ceptually midway between the Klimontovich equation and fluid models. It de- scribes the time evolution of the distribution function f s (x, v, t) in so-called phase space - space with spatial x and velocity v dimensions [3].

∂f

∂t + v · ∇ x f + q

m [E + v × B] · ∇ v f = 0. (2) The computational method of this thesis focuses on solving the Vlasov’s equation in combination with electromagnetic description.

Vlasov-Maxwell System of Equations forms a self-consistent plasma de-

scription in 3D space by combining the Vlasov equation (2) and Maxwell’s

equations (3):

(8)

 

 

 

 

 

 

 

 

∇ · E = ρ

 0

,

∇ × E = − ∂B

∂t ,

∇ · B = 0,

∇ × B = µ 0

 j +  0

∂E

∂t

 .

(3)

The coupling between the Vlasov equation (2) and the Maxwell system of equations (3) is achieved through fluid quantities - charge density ρ and current density j:

ρ(x, t) = X

s

q s

Z

f (x, v, t) dv, (4)

j(x, t) = X

s

q s

Z

vf (x, v, t) dv. (5)

Electrostatic limit is used in the absence of the magnetic field and Maxwell’s equations (3) reduce to the Gauss’ law:

∂E

∂x = X

s

q s n s . (6)

This formulation that consists of the Vlasov equation (2) and Gauss law (6) is called Vlasov-Poisson system of equations. In this thesis electrostatic formu- lation is used for 1D1V simulations and fully electromagnetic Vlasov-Maxwell for 2D3V simulations.

1.1 Numerical Methods for Plasma Physics

The three most common computational approaches to solve the Vlasov-Maxwell System of Equations are grid based methods (Eulerian approach), particle-based methods and spectral methods.

The most common technique is the particle-based Particle-in-Cell (PIC) method [4, 5, 6]. The PIC method statistically represents the distribution func- tion by means of computational super-particles. Because of the limited number of particles used in simulations, the PIC method is subject to numerical noise that effects the accuracy of the simulations. A large number of particles must be used to have a good representation of the particle distribution function and reduce noise that comes from the discrete nature of particles [7].

Eulerian grid-based methods discretize phase space on a multidimensional grid. For solving a typical kinetic problem - electromagnetic wave interaction with electrons very fine 6D spatial and velocity space discretization and short time steps must be used in order to resolve the smallest time and spatial scales in the system - electron gyroradius, gyrofrequency or plasma frequency.

Spectral methods are alternative to co-existing particle and grid based meth- ods and represent distribution function with continuous basis functions in veloc- ity and spatial space. Spectral methods like Hermite-Fourier (HF) removes the problem of statistical noise, allowing for kinetic simulations with high accuracy.

It can be used in plasma regimes that are not practical to solve with grid-based

(9)

and PIC methods. Spectral methods for plasma physics are complex and less robust than particle and grid-based methods, but are of very high interest for simulations of weak turbulence and wave-wave interaction.

1.2 Scope of Thesis

This thesis focuses on the Hermite-Fourier spectral method which discretizes the velocity space in Hermite basis and the spatial space in Fourier basis.

The goal of this thesis is to improve robustness and convergence of the method, parallelize existing code and study its performance.

This thesis is organized as follows. In the second section, the HF method is

described and the governing equations are given. Derivations start from simple

1D1V case and end with a multi-dimensional 2D3V description. The third

section presents novel improvements of the HF method such as dynamically

adaptive change of Hermite basis and use of high-end scientific libraries. The

fourth section focuses on parallelization of the code on shared memory systems

and performance tuning. Finally, in the fifth section a summary is given and

future work is discussed.

(10)

2 Literature Survey and Theoretical Background

This section discuses the development of the Hermite-Fourier (HF) spectral method from its initial proposal to its recent extension to multidimensional case. At the beginning a 1D electrostatic method is presented, it illustrates the procedure of transformation from physical space to spectral one, solution of the system of partial differential equations (PDE) and transformation back to real space. Later in this section the full extension to multidimensional case is shown.

2.1 The Hermite-Fourier Spectral Method

The use of Hermite functions to solve the Vlasov equation was proposed for the first time by Grad [8] in 1949. It was only in 1998 that an explicit numerical solution of the Vlasov-Poisson system with the HF method was obtained by Schumer [9] in one dimension and in electrostatic limit. Since then the method has been improved in terms of formulation, numerical discretization and recently its convergence and conservation properties have been studied. For example the Crank-Nicolson scheme to discretize the problem in time has been presented in [10], it was shown that this scheme leads to an exact conservation of the total mass, momentum and energy. Conservation of all three physical quantities is challenging or even impossible to achieve with many other simulation methods.

2.2 Velocity Expansion in Hermite Basis

The main idea of the HF spectral method is to transform distribution function and the Vlasov-Poisson system of equations to the Hermite-Fourier spectral space. Firstly, we start by normalizing the system to non-dimensional quanti- ties, see Nomenclature section. For simplicity we remove “ˆ” from normalized quantities, then the Vlasov equation for the species s in presence of the electric field E(x, t) is

∂f s

∂t + v ∂f s

∂x + q s m s

E ∂f s

∂v = 0. (7)

In the electrostatic limit, the Maxwell equations reduce to Gauss’ law:

∂E

∂x = X

s

q s n s , (8)

where n s = R f s (x, v, t) dv is density of the species s. Equation (7) and (8) form Vlasov-Poisson equation system and self-consistent representation of 1D electrostatic plasma. This system is transformed to Hermite basis.

There are two kinds of Hermite basis function - symmetrically and asymmet- rically weighted. This thesis follows [10, 11, 12] and uses the Asymmetrically- Weighted (AW) Hermite functions to describe the velocity dependence of the distribution function. The AW Hermite basis functions are defined as:

Ψ n (ξ) = (π2 n n!) −1/2 H n (ξ)e −ξ

2

,

Ψ n (ξ) = (2 n n!) −1/2 H n (ξ), (9)

where H n is the Hermite polynomial of the n-th order [13]. The first six basis

functions are shown in Figure 1. It is important to note in Figure 1 that the 0-th

(11)

order AW Hermite polynomial is a Maxwellian function. For this reason, just one AW Hermite polynomial is sufficient to describe the distribution function of a Maxwellian plasma.

Figure 1: First six Hermite basis functions Ψ n (ξ).

Normalized distribution function defined in Hermite spectral space as f ˆ s (x, v, t) =

N

H

−1

X

n=0

C n s (x, t)Ψ n (ξ s ), (10)

where C n s are expansion coefficients of Hermite mode n and species s. In order to improve convergence of the Hermite basis, it can be shifted and scaled

ξ s = (v − u s )

α s , (11)

where u s and α s are normalized (non-dimensional) shift and scaling parameters.

Three useful relations of the Hermite basis functions [14] are the orthonor- mality condition

Z ∞

−∞

Ψ n (ξ)Ψ m (ξ) dξ = δ n,m , (12) the recursion relation

n (ξ) =

r n + 1

2 Ψ n+1 (ξ) + r n

2 Ψ n−1 (ξ) (13)

and the derivative relation dΨ n (ξ)

dξ = − p

2(n + 1)Ψ n+1 (ξ). (14)

For scaled and shifted basis (11) recursion and derivative relations (13) and (14) become

ˆ

n (ξ) = α s

r n + 1

2 Ψ n+1 (ξ) + α s r n

2 Ψ n−1 (ξ) + u s Ψ n (ξ), (15) dΨ n (ξ)

dˆ v = − 1 α s

p 2(n + 1)Ψ n+1 (ξ). (16)

(12)

Now follows actual transformation of the Vlasov-Poisson system to the Hermite- Fourier basis. The distribution function (10) is inserted in (7), multiplied by Ψ m and integrated over dξ. Summation series are finite, hence one can interchange summation and integration signs:

N

H

−1

X

n=0

∂C n s (x, t)

∂t

Z ∞

−∞

Ψ n (ξ)Ψ m (ξ) dξ +

N

H

−1

X

n=0

∂C n s (x, t)

∂x

Z ∞

−∞

n (ξ)Ψ m (ξ)dξ

+ q s

m s

E

N

H

−1

X

n=0

C n s (x, t) Z ∞

−∞

∂Ψ n (ξ)

∂v Ψ m (ξ)dξ = 0. (17) Now we apply recursion (15), derivative (16) and orthonormality (12) rela- tions:

N

H

−1

X

n=0

∂C n s (x, t)

∂t δ n,m

+

N

H

−1

X

n=0

∂C n s (x, t)

∂x α s

r n + 1

2 δ n+1,m + α s r n

2 δ n−1,m + u s δ n,m

!

+ q s m s

E

N

H

−1

X

n=0

C n s (x, t)



− 1 α

p 2(n + 1)δ n+1,m



= 0. (18)

From Kronecker-delta property P

n

a n δ n,m = a m follows

∂C m s

∂t +α s

r m + 1 2

∂C m+1 s

∂x + r m 2

∂C m−1 s

∂x + u s

α s

∂C m s

∂x

!

− q s

m s

√ 2m

α s E C m−1 s = 0.

(19) Similarly, the Gauss law is transformed to the Hermite basis. By inserting (10) and (9) in (8) one gets

∂E

∂x = X

s

α s q s N

H

−1

X

n=0

C n s (x, t)(π2 n n!) −1/2 Z ∞

−∞

H n (ξ)e −ξ

2

dξ, (20)

where dv = α s dξ s follows from (11). The integral R ∞

−∞ H n (ξ)e −ξ

2

ξ for n = 0 is non-zero, for n > 0 it is zero, consequently,

∂E

∂x = X

s

q s α s C 0 s (x). (21)

2.3 Spatial Expansion in Fourier Basis

Fourier basis is the most popular choice for expansion of the distribution func- tion in space [15]. Waves in plasma can be represented as a series of periodic and continuous Fourier basis functions Φ k (x):

Φ k (x) = exp  2πikx L



. (22)

(13)

Distribution function (10) which was transformed to Hermite basis now is also expanded in Fourier spatial basis:

f s (x, v, t) =

N

H

−1

X

m=0 N

X

k=−N

C m,k s (t)Ψ msk (x). (23)

Expansion coefficients

C m s (x) =

N

X

k=−N

C m,k s e

2πikxL

, (24)

E(x) =

N

X

k=−N

E k e

2πikxL

(25)

are inserted in (19)

N

X

k=−N

∂C m,k s

∂t e

2πikxL

+ 2πikα s

L

r m + 1 2

N

X

k=−N

C m+1,k s + r m 2

N

X

k=−N

C m−1,k s + u s

α s N

X

k=−N

C m,k s

! e

2πikxL

− q s

m s

√ 2m α s

N

X

k=−N

E k e

2πikxL

N

X

k=−N

C m−1,k s e

2πikxL

= 0. (26)

Useful properties of the Fourier basis are orthogonality condition Z L

0

exp  2πikx L

 exp



− 2πinx L



= Lδ k,n (27)

and Cauchi product for finite series

N

X

k=−N

a k

! N X

k=−N

b k

!

=

N

X

k=−N N

X

j=−N

a j b k−j . (28)

Cauchi product is multiplication of two finite series and the result is discrete convolution of these series.

Equation (26) is multiplied by Φ −n (x) (see (22)) and Cauchi product is used (28), finally, expression is integrated over 0 ≤ x ≤ L:

Z L 0

N

X

k=−N

dC m,k s

dt e

2πikxL

e

2πinxL

dx

+ Z L

0

α s

r m + 1 2

N

X

k=−N

C m+1,k s + r m 2

N

X

k=−N

C m−1,k s + u s α s

N

X

k=−N

C m,k s

!

e

2πikxL

2πik

L e

2πinxL

dx

− Z L

0

q s m s

√ 2m α s

N

X

k=−N k

X

j=−N

C m−1,j e

2πijxL

E k−j e

2πi(k−j)xL

e

2πinxL

dx = 0. (29)

(14)

After applying orthogonality condition (27) and property of Dirac delta a set of ordinary differential equations for expansion coefficients C n,k s (t) is obtained:

dC m,n s dt + α s

2πin L

r m + 1

2 C m+1,n s + r m

2 C m−1,n s + u s

α s C m,n s

!

− q s

m s

√ 2m α s

N

X

j=−N

E n−jC

m−1,j

= 0. (30)

Similarly, Gauss law is transformed by substituting (25) in (21)

N

X

k=N

E k e

2πikxL

2πik

L = X

s

α s q s N

X

k=N

C 0,k s e

2πikxL

. (31)

The result is multiplied by Φ −n , integrated over dx and after applying the orthogonality condition (27) we get relation between expansion coefficients of the distribution function and electric field

E k

2πik

L = X

s

q s α s C 0,k s , (32)

where the result for k = 0 in not defined. Relation is rewritten in simpler form:

E k =

0, if k = 0

− iL 2πk

P

s

q s α s C 0,k s , otherwise. (33)

2.4 Filamentation Problem and Collisional Operator

A common problem for all spectral solutions of the Vlasov equation is filamen- tation of phase space [16] - the development of increasingly small phase space structures that cannot be resolved within the number of spectral modes in use.

The filamentation has its typical time scale called “recurrence time”. After this time the system returns to its initial state. The recurrence time is proportional to the square root of the number of the Hermite modes used in the simulation.

For this reason, increasing the number of Hermite modes increases the recur- rence time but it does not eliminate the filamentation problem. One of the numerical techniques to address this problem is to add a collisional term to the right hand of (30). Such technique is implemented in [17, 11]. For this purpose, the Lenard-Bernstein collision operator C LB [17] is used:

C LB [ ˆ f ] = ˆ ν ∂

∂ ˆ v ˆ v ˆ f + ∂ ˆ f

∂ ˆ v

!

. (34)

Collisional operator the is modified according to [11]¿

C[C m,k ] = ν m(m − 1)(m − 2)

(N H − 1)(N H − 2)(N H − 3) C m,k , (35)

where ν is collisionality parameter, m and k are Hermite and Fourier modes

respectively. Equation (30) with added collisional term is

(15)

dC m,n s

dt + α s 2πin L

r m + 1

2 C m+1,n s + r m

2 C m−1,n s + u s α s

C m,n s

!

− q s m s

√ 2m α s

N

X

j=−N

E n−j C m−1,j s + ν m(m − 1)(m − 2)

(N H − 1)(N H − 2)(N H − 3) C m,n s = 0.

(36)

2.5 Implicit Time Discretization

Recent works propose semi-implicit and fully implicit time discretizations of the HF method [10, 11]. It is shown that fully implicit time discretization leads to an exact conservation of the total mass, momentum and energy.

Applying the Crank-Nicolson scheme to equation (36) we get C m,n t+1 − C m,n t

∆t +α s

2πin L

r m + 1 2

(C m+1,n t+1 + C m+1,n t )

2 + r m

2

(C m−1,n t+1 + C m−1,n t )

2 + u s

α s

(C m,n t+1 + C m,n t ) 2

!

− q s

m s

√ 2m α s

N

X

j=−N

(E n−j t+1 + E n−j t ) 2

(C m−1,j t+1 + C m−1,j t ) 2

+ ν m(m − 1)(m − 2) (N H − 1)(N H − 2)(N H − 3)

"

C m,n t+1 + C m,n t 2

#

= 0. (37)

2.6 Transformation to Real Space

Definition of distribution function (23) in the Hermite-Fourier spectral space is used to reconstruct solution in the real space. By integrating (23) over velocity space and neglecting spatial expansion in the Fourier basis, we obtain

n(x) = Z

f (x, v) dv =

Z N

H

−1 X

n=0

C n (x)Ψ n (v) dv = C 0 (x). (38)

Similarly, the fluid (or bulk) velocity is defined as v avg = n 1 R vf dv and we can express it in terms of Hermite coefficients as

v avg = R vf dv

n =

Z v

N

H

−1

X

n=0

C n (x)Ψ n (v) dv = C 1 (x)

√ 2C 0 (x) . (39)

The thermal velocity is defined as v th = r R (v − v avg ) 2 f dv

n and transformed into the Hermite space:

v th =

r R (v − v avg ) 2 f dv

n = α(x)

√ 2 s

1 + √ 2 C 2 (x)

C 0 (x) − C 1 (x) 2

C 0 (x) 2 . (40)

(16)

The pressure stress tensor P in the laboratory frame is

P = Z

mv 2 f dv = Z

mv 2

N

H

−1

X

n=0

C n (x)Ψ n (ξ(v)) dv = 1 2 m 

C 0 (x) + √

2C 2 (x)  . (41) The kinetic energy for each species is written as

W kin = 1 2

Z L 0

Z ∞

−∞

mv 2 f (x, v, t)dxdv = Lmα(2u 2 + α 2 )

4 C 0,0 + Lmuα 2

√ 2 C 1,0 + Lmα 3 2 √

2 C 2,0 . (42) The electric field energy is

W el = 1 2

Z L 0

E 2 dx = L 2

N

X

k=−N

|E| 2 . (43)

2.7 Extension to the Multi-Dimensional Case

Extension of the HF method to multiple dimensions (3D3V) and to the fully electromagnetic (E,B) case is given by Delzanno [12]. In this section notation and general equations used in the paper are given.

The Vlasov-Maxwell system is normalized to non-dimensional quantities, see Nomenclature section. For simplicity we remove “ˆ” from normalized quantities, then the Vlasov-Maxwell system for the species s in presence of the electric E(x, t) and magnetic B(x, t) fields is

∂f s

∂t + v · ∇ x f s + q s

e ω cs

ω pe

(E + v × B) · ∇ v f s = 0, (44)

∂B

∂t = −∇ × E, (45)

∂E

∂t = ∇ × B − ω pe ω ce

X

s

 q s

e

+∞

Z

−∞

vf s dv

 , (46)

∇ · E = ω pe

ω ce X

s

 q s

e

+∞

Z

−∞

vf s dv

 , (47)

∇ · B = 0. (48)

Equations (9) and (22) that were used for 1D1V discretization are also used for 2D3V electromagnetic case. Distribution function in 2D3V Hermite-Fourier spectral space is

f s (x, v, t) =

N

n

−1

X

n=0 N

m

−1

X

m=0 N

p

−1

X

p=0

C n,m,p k

x

,k

y

,s (t) Ψ n (ξ x s )Ψ m (ξ s y )Ψ p (ξ z s )

N

x

X

x=−N

x

N

y

X

y=−N

y

exp



2πi  k x x L x

+ k y y L y



. (49)

(17)

To transform the Vlasov equation to Hermite space, the distribution function (49) is inserted in the Vlasov’s equation (44), multiplied by orthogonal Hermite basis Ψ nx smy spz s ), integrated over dξ x sy sz s and the orthogonality condition (12) is applied. Then to transform to the Fourier space, the result is multiplied by exp



−2πi  k 0 x x L x + k y 0 y

L y



and integrated over spatial domain.

The final system of equations for expansion coefficients in spatial and velocity spaces is

dC n,m,p k

x

,k

y

,s

dt = − 2πik x L x

α s x

r n + 1

2 C n+1,m,p k

x

,k

y

,s + r n

2 C n−1,m,p k

x

,k

y

,s + s u s x

α x s C n,m,p k

x

,k

y

,s

!

− 2πik y L y

α y s

r m + 1

2 C n,m+1,p k

x

,k

y

,s + r m

2 C n,m−1,p k

x

,k

y

,s + s u s y

α s y C n,m,p k

x

,k

y

,s

!

+ q s

e ω pe

ω ce

√ 2n

α s x [E x ∗ C n−1,m,p ] k

x

,k

y

+

√ 2m

α s y [E y ∗ C n,m−1,p ] k

x

,k

y

+

√ 2p

α s z [E z ∗ C n,m,p−1 ] k

x

,k

y

!

+ q s e

ω pe ω ce

√ mp  α s z α s y − α s y

α s z



[B x ∗C n,m−1,p−1 ] k

x

,k

y

+ p

m(p + 1) α s z

α s y [B x ∗C n,m−1,p+1 ] k

x

,k

y

− p

(m + 1)p α s y

α s z [B x ∗C n,m+1,p−1 ] k

x

,k

y

+ √ 2m u s z

α s y [B x ∗C n,m−1,p ] k

x

,k

y

− p 2p u s y

α s z [B x ∗C n,m,p−1 ] k

x

,k

y



+ q s

e ω pe

ω ce

√ np  α s x α s z − α s z

α s x



[B y ∗C n−1,m,p−1 ] k

x

,k

y

+ p

p(n + 1) α s x

α z s [B y ∗C n+1,m,p−1 ] k

x

,k

y

− p

(p + 1)n α s z

α s x [B y ∗C n−1,m,p+1 ] k

x

,k

y

+ p 2p u s x

α s z [B y ∗C n,m,p−1 ] k

x

,k

y

− √ 2n u s z

α s x [B y ∗C n−1,m,p ] k

x

,k

y



+ q s e

ω pe ω ce

 √ nm

 α s y α s x − α x s

α s y



[B z ∗C n−1,m−1,p ] k

x

,k

y

+ p

n(m + 1) α s y

α s x [B z ∗C n−1,m+1,p ] k

x

,k

y

− p

(n + 1)m α s x

α s y [B z ∗C n+1,m−1,p ] k

x

,k

y

+ √ 2n u s y

α s x [B z ∗C n−1,m,p ] k

x

,k

y

− √ 2m u s x

α y s [B z ∗C n,m−1,p ] k

x

,k

y

 . (50) Similarly, Maxwell’s equations transformed to the Fourier space result in a

set of six equations are dE x k

x

,k

y

dt = 2πi k y L y

B z k

x

,k

y

− ω pe ω ce

X

s

q s e



α s x α s y α z s  α s x

√ 2 C 1,0,0 k

x

,k

y

,s + u s x C 0,0,0 k

x

,k

y

,s



, (51) dE y k

x

,k

y

dt = −2πi k x L x

B z k

x

,k

y

− ω pe ω ce

X

s

q s e

 α s x α y s α s z

 α s y

√ 2 C 0,1,0 k

x

,k

y

,s + u s y C 0,0,0 k

x

,k

y

,s



, (52) dE z k

x

,k

y

dt = 2πi  k x L x

B y k

x

,k

y

− k y

L y

B k x

x

,k

y



− ω pe

ω ce

X

s

q s

e



α s x α s y α s z  α s z

2 C 0,0,1 k

x

,k

y

,s + u s z C 0,0,0 k

x

,k

y

,s



, (53)

dB x k

x

,k

y

dt = −2πi k y

L y E z k

x

,k

y

, (54)

(18)

dB k y

x

,k

y

dt = 2πi k x

L x

E z k

x

,k

y

, (55)

dB k z

x

,k

y

dt = −2πi  k x

L x E y k

x

,k

y

− k y

L y E k x

x

,k

y



. (56)

(19)

3 Development of the Hermite-Fourier Method

In this section improvements of the HF spectral method are presented. Firstly, a technique of dynamically adaptive number of Hermite modes is discussed, the purpose of this technique is to increase efficiency of large-scale simulations. Sec- ondly, transformation method that increases convergence and robustness of the method is discussed. Finally, implementation of non-linear solver and precon- ditioners from PETSc library is presented.

3.1 Dynamically Adaptive Number of Fluid Moments

In Section 2.6 it was shown how fluid quantities like density, bulk velocity and pressure are calculated. It can be verified that each additional fluid moment will require expansion coefficients of higher Hermite mode. In [18] it is shown that solving the problem in HF space is equivalent to solving fluid equations.

For example, coefficients C 0 and C 1 obtained from equations (38) and (39) in terms of fluid quantities

C 0 = n ⇒ C 1 = √

2v avg n (57)

and inserting in the first equation (m=0) of the system (19) which is dC 0,k s

dt + α s

2πik L

r 1

2 C 1,k s + u s

α s

C 0,k s

!

= 0 (58)

continuity equation is obtained

∂n

∂t + ∂

∂x (v avg n) = 0. (59)

From equations (41) coefficient C 2 is C 2 = − mn − 2P

√ 2m . (60)

By inserting (57) and (60) in the second equation (m=1) of the system (19) which is

∂C 1 s

∂t + α s

∂C 2 s

∂x + r 1

2

∂C 0 s

∂x + u s

α s

∂C 1 s

∂x

!

− q s

m s

√ 2 α s

E C 0 s = 0 (61)

equation for conservation of momentum is obtained for u = 0 and α = 1:

∂t (mn v avg ) + ∂P

∂x − qnE = 0. (62)

Similarly, by substituting higher expansion coefficients in the equation sys-

tem (19) conservation equations are obtained - energy and higher moment equa-

tions.

(20)

3.1.1 Dynamically Adaptive Change of the Number of Hermite Modes Plasmas in presence of collisions tend to relax to Maxwellian velocity distribu- tion. For Maxwellian plasmas fluid description can be used. Only four modes are required for fluid description of plasma to conserve mass, momentum and energy. Collisionless plasmas on the other hand may have very complex velocity distribution and kinetic description of plasma is required.

Adding or removing Hermite modes provides smooth transition between fluid and kinetic descriptions. In [18] it is shown that the number of Hermite modes can be changed during a simulation. Hermite polynomials are added/removed for species s according to this rule:

- add one Hermite mode if max(|C N s

H

−1,k |) > C 0,0 s H tol + - remove one Hermite mode if max(|C N s

H

−1,k |) < C 0,0 s H tol , (63) where C 0,0 is the 0-th order expansion coefficient which is proportional to the average density, H tol + and H tol are the threshold values, and max represents the maximum value among all Fourier modes.

3.1.2 Results

In [18] two tests are carried to demonstrate the method. Initial conditions for Landau damping with ions and two-stream instability simulations are given in Table 1.

Table 1: Initial parameters for the Landau damping with ions and the two stream instability simulations.

Param. Land. damp. w. ions Two str. inst. Description

L 2π 2π Simulation domain length

T 70 70 End time

∆t 0.05 0.05 Time step

N F 25 25 Number of Fourier modes

v th 1 & 1/ √

1836 1 Thermal velocity

ν 10 50 Collisional parameter

u avg 0 ±1 Average velocity

n 1 & 1 0.5 + 0.5 Average density

qom −1 & 1/1836 −1 Charge over mass ratio

 0.001 0.001 Initial perturbation

H tol + 10 −12 10 −12 Upper threshold value

H tol 10 −14 10 −14 Lower threshold value

Landau Damping With Ions. Landau damping is an effect of longitudinal

wave damping due to collisionless interaction between charged particles. For

electrons these longitudinal waves are called Langmuir waves and for ions these

waves are called ion acoustic waves. Usually ions are neglected in Landau damp-

ing simulations because of high mass. Ion thermal velocity is much lower than

electron thermal velocity therefore the damping rate for an ion acoustic wave is

much smaller than for electrons.

(21)

Parameters used in the simulation are given in Table 1. Results are shown in Figure 2 (left). The black solid line shows the first mode of the electric field E 1 . Theoretical damping rate is given by black dashed line. E 1 shows two electric field modes present in plasma - Langmuir and ion acoustic. At the beginning of simulation only electrons are perturbed. Electrons interact with ions and start an ion acoustic wave with a much smaller amplitude. Since the damping rate of the Langmuir wave is larger than the damping rate of the ion acoustic wave, the second is present after the Langmuir wave damps out. The number of Hermite modes for ions (red solid line) increases at the beginning of the simulation and saturates at only 7 modes, while for electrons (red dotted line) the number of modes for electrons increases initially and reaches the peak value of 139 modes. This clearly shows the benefit of using an adaptive number of modes for multiscale simulations.

Two Stream Instability can be initiated by two oppositely-directed beams.

While the instability is growing, beam kinetic energy is transformed to potential energy.

Parameters used in the simulation are given in Table 1. Results are shown in Figure 2 (right). The black solid line shows the first mode of the electric field E 1 . Theoretical damping rate is given by the black dashed line. The number of electron modes (red solid line) grows slowly during the linear phase (at times t = 0 − 40 ω pe ), but increases and saturates up to 195 modes in nonlinear phase (at times t = 40 − 70 ω pe ).

0 10 20 30 40 50 60 70

10−9 10−8 10−7 10−6 10−5 10−4 10−3

t (1/ω

p

)

|E

1

|

0 10 20 30 40 50 60 70

0 20 40 60 80 100 120 140

t (1/ω

p

)

N

H Ions Electrons

0 10 20 30 40 50 60 70

10−7 10−6 10−5 10−4 10−3 10−2 10−1 100 101

t (1/ω

p

)

|E

1

|

0 10 20 30 40 50 60 70

0 20 40 60 80 100 120 140 160 180 200

t (1/ω

p

)

N

H

Figure 2: Landau damping (left) and two-stream instability (right) tests. The black solid line represents the amplitude of the first spectral mode of the electric field E 1 (left black axis). Black dashed lines represent damping and growth rates from the linear theory. Red lines represent the number of Hermite modes for different species. Pictures taken from [18].

Performance. Adaptive number of Hermite modes can reduce the number of equations that need to be solved and save computational time. In the paper [18] savings in terms of Hermite modes are defined as following:

% N H = P

s

(N H ref s − N H adapt s ) P

s

N H ref s · 100, (64)

(22)

where N H ref s and N H adapt s are the number of Hermite modes for the refer- ence and adaptive simulations for particle species s. To compare improvement two simulations are run, first adaptive and second simulation with fixed num- ber of modes. For the reference simulation the number of Hermite modes is set to the maximum number obtained for the adaptive simulation: N H ref s = max(N H adapt s ). The parameter is computed at each simulation cycle and rep- resents theoretical time savings from the adaptive approach. Similarly, the computational savings are computed as:

% CP U = t ref − t adapt

t ref

· 100, (65)

where t ref and t adapt are computation time for reference and adaptive simula- tions.

Efficiency of the technique is shown in Figure 3. Reduced number of modes requires less linear equations to solve. Results show that CPU time per time step strongly correlates with the number of Hermite modes.

Landau damping simulation Figure 3 (left) with the adaptive number of Hermite mode takes 46.7 s while the reference simulation that uses fixed number of Hermite modes takes 56.3 s. The adaptive technique leads to a 21% overall time reduction.

For the two-stream instability test Figure 3 (right) the CPU saving factor matches the Hermite mode saving factor. At the beginning of the simulation a reduced number of modes allows the simulation to run about 10 times faster than the reference simulation. The total computation time for adaptive simulation is 58 s versus 110 s for reference simulation. This leads to 47% overall time reduction compared to the reference simulation time.

0 10 20 30 40 50 60 70

0 10 20 30 40 50 60 70 80 90 100

t (1/ ω

p

)

%

CPU time saved NH saved

0 10 20 30 40 50 60 70

0 10 20 30 40 50 60 70 80 90 100

t (1/ ω

p

)

%

CPU time saved NH saved

Figure 3: Savings in terms of time and Hermite modes for technique of adaptive number of Hermite modes. Landau damping (left) and two-stream instability (right) tests. Pictures taken from [18].

3.2 Dynamically Adaptive Change of Hermite Basis

In this section the method for changing velocity scaling and shifting parameters

of the Hermite basis during a simulation is described. The method is based on

the physical interpretation of velocity scaling and shifting parameters that are

proportional to average and thermal velocities, respectively.

(23)

Plasma velocity distribution can change considerably during simulation, for example, perturbed isotropic plasma tends to relax to the equilibrium state which is Maxwellian distribution and results in plasma heating. Plasma beams in the presence of background plasma can lose energy and slow down. In cases where plasma thermal and average velocities change significantly, the use of right parameters can result in a substantial reduction of Hermite modes and computational savings.

3.2.1 Basis Transformation Method

In order to change the Hermite basis, expansion coefficients must be transformed from the old basis to the new one. The distribution function in real space is written as two different expansions in Hermite basis (10) and Fourier basis (24):

f (x, v, t) =

N

H

−1

X

n=0 N

X

l=−N

C n,l old (t)Ψ n (ξ old )e

2πilxL

=

N

H

−1

X

m=0 N

X

k=−N

C m,k new (x, t)Ψ m (ξ new )e

2πikxL

. (66) Both sides are multiplied by Ψ mnew )e

2πik0 xL

that are orthogonal functions to the new basis. After integrating over dξ new dx we get expansion coefficients in the new basis:

L C m,k new = Z ∞

−∞

N

H

−1

X

n=0 N

X

l=−N

C n,l old (t)Ψ n (ξ old )e

2πilxL

Ψ m (ξ new )e

2πik0 xL

dξ new dx.

(67) Parameters u and α are averaged in space and chosen independently from coordinate x, which means that ξ is independent from x, and the same transfor- mation procedure can be applied on all Fourier modes. Orthogonality of Fourier series (27) is used to simplify the previous equation

C m,k new = Z ∞

−∞

N

H

−1

X

n=0

C n,k old (t)Ψ noldmnew ) dξ new . (68) Summation and integration can be swapped and transformation of expansion coefficients can be written in the form

C n,k new = T n,m C m,k old . (69) To obtain the transformation matrix, we solve the integral

T n,m = Z ∞

−∞

Ψ n (ξ old )Ψ m (ξ new ) dξ new . (70)

The integral from (70) can be simplified by assuming that basis transformation can be applied independently for shift and scaling parameters. Transformation matrix T n,m can be split into two independent transformations - one for scaling and another for shifting:

T n,m = T n,h α T h,m u . (71)

(24)

Transformation for shift parameter is obtained by setting α = const and solving (70). This integral was solved by means of symbolic software and the result is a general expression for constructing the transformation matrix for shift parameter:

T h,m u = s

m! 2 m−h h! ((m − h)!) 2

 u old − u new α

 m−h

, (72)

where symbol ”!” represents factorial.

Transformation for scaling parameter is obtained in a similar way. The drift velocity is set to u = const and equation (70) is solved. The general expression for constructing the transformation matrix for scaling parameter is:

T n,h α =

 

 

0, if h + n odd

s

h!! (h − 1)!!

n! ((h − n)!!) 2 α n+1 old α h+1 new

2 old − α 2 new ) (h−n)/2 , if h + n even (73)

where symbol ”!!” represents double factorial.

The choice of shifting and scaling parameters has been studied by many authors. How different scaling parameters affect the convergence of solution was studied by [9], while the more mathematical technique of finding the optimal pa- rameters in orthogonal systems described by [19] and for Hermite basis applied by [20].

In this thesis a qualitative rather than quantitative understanding of opti- mal choice of parameters is used. Transformations (72) and (73) are used to study implications of not optimally chosen basis parameters. The same ini- tial distribution function which is simple Maxwellian is transformed in different basis.

In Figure 4 basis shifting for a simple Maxwellian-like distribution is demon- strated. Initially the distribution is expressed with one mode C 0 = 1 in the basis with α = 1 and u old = 0. Then, by using transformation matrix (72) the initial distribution is transformed in the new basis where u new = 1 (left panel) and u new = 2 (right panel). The distribution function from expansion in the new basis converge for large number of Hermite modes.

In Figure 5 basis scaling for a simple Maxwellian-like distribution is demon-

strated. Initially the distribution is expressed with one mode C 0 = 1 in the

basis with α old = 1 and u = 0. Then, by using transformation matrix (73) the

initial distribution is transformed in the new basis where α new = 2 (left panel)

and u new = 0.65 (right panel). It was found that the distribution function can

be expanded in slightly scaled basis with sufficient number of Hermite modes,

although scaling to a basis with much smaller α in addition to a large num-

ber of modes may result in a loss of positivity and an overall unsatisfactory

representation of the initial distribution.

(25)

original 3 modes 7 modes 13 modes

-2 -1 1 2 v

0.1 0.2 0.3 0.4 0.5 0.6 0.7 f(v)

original 5 modes 13 modes 23 modes

-2 -1 1 2 v

-0.5 0.5 1.0 f(v)

Figure 4: Basis shift example. Blue solid line (original) represents the dis- tribution function for the expansion coefficient C 0 = 1 in basis with α = 1 and u old = 0. Dashed coloured lines represent distribution functions for speci- fied number of Hermite modes in shifted basis with u new = 1 (left panel) and u new = 2 (right panel).

original 3 modes 7 modes 13 modes

-2 -1 1 2 v

0.1 0.2 0.3 0.4 0.5 f(v)

original 5 modes 13 modes 23 modes

-2 -1 1 2 v

-0.5 0.5 1.0 f(v)

Figure 5: Basis scaling example. Blue solid line (original) represents the dis- tribution function for the expansion coefficient C 0 = 1 in basis with α old = 1 and u = 0. Dashed coloured lines represent distribution functions for speci- fied number of Hermite modes in scaled basis with α new = 2 (left panel) and α new = 0.65 (right panel).

Distribution function errors from improper choice of basis are two type - numerical errors and errors from finite expansion in Hermite basis. Fig- ures 4 and 5 show that errors from finite expansion in the Hermite basis are less intuitive, hence similar basis transformations are performed and shown in Figures 4 and 5 for different scaling and shift parameters.

The error for two distribution functions is calculated as

 = v u u u t

v

max

Z

v

min

[f orig (v) − f transf (v)] 2 dv, (74)

where f orig and f transf are original Maxwellian distribution and distribution obtained by basis transformation. The integral is computed numerically and results are shown in Figure 6. For basis shifting (left) error is larger for large shift parameters and small number of Hermite modes. Results for 75 and 100 modes intersect which is believed to be a result of limited numerical precision.

For basis scaling (right) note that the error behaves differently for parameters

smaller and larger than 1. For larger coefficients it is expected that the error

(26)

asymptotically approach some limit, while for scaling coefficients smaller than 1 error is not limited and for α / 0.7 more Hermite modes introduce larger error due to oscillations shown in Figure 5 (left).

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

10−30 10−25 10−20 10−15 10−10 10−5 100 105

u

ε

25 50 75 100

0.6 0.8 1 1.2 1.4 1.6 1.8 2

10−30 10−20 10−10 100 1010 1020 1030

α

ε

25 50 75 100

Figure 6: Distribution function error  as a function of shift parameter u (left) and scaling parameter (right) for different number of Hermite modes (lines).

Implementation of the method. Physical interpretation of velocity scal- ing and shifting parameters means that both are proportional to average and thermal velocities, respectively. The optimal shift parameter is defined as the average velocity (39):

u new = v avg = u old + α

√ 2 avg x  C 1 (x) C 0 (x)



. (75)

The optimal scaling parameter is defined as the average thermal velocity (40)

α new = √

2v th = α old avg x

"s 1 + √

2 C 2 (x)

C 0 (x) − C 1 (x) 2 C 0 (x) 2

#

, (76)

where factor √

2 before v th follows from Maxwellian velocity distribution and definition of α.

Quantities C 0 (x), C 1 (x) and C 2 (x) are expansion coefficients in Hermite basis and in real spacial space, these coefficients are calculated by (24). One should be aware that C 0 (x) is proportional to the particle density and numerical exception (invalid floating point operation) can occur if this quantity is small.

Change of basis requires multiplication of large matrices, therefore we in- troduce tolerances for parameters to update basis when it is needed. We do transformation of basis for the shift parameter when

u tol < |u old − u new |, (77) where u tol is absolute torelance for shift parameter. Similarly, we do transfor- mation of scaling parameter when

α tol < |α old − α new | α old

, (78)

where α tol is relative tolerance for scaling parameter. The right choice of these

tolerances is problem and implementation specific. Each transformation of basis

(27)

requires generation of the transformation matrix, transformation of expansion coefficients and update of Hermite basis functions (9), hence it is justified to not to set both to zero. After performing many tests it was found that small and large tolerances may result in numerical underflow and overflow, hence it is advisable to keep booth tolerances bounded.

3.2.2 Results

In this section results from two simulations are present: the non-linear Landau damping and the two-stream instability. In both simulation tests a wave is triggered by perturbing the electron distribution function

f e (x, v) = n 0 [1 +  cos(2πx/L)]

√ πα e

exp

"

−  v − u e α e

 2 #

, (79)

where  is the amplitude of the initial perturbation. The initial conditions for both simulations are shown in Table 2.

Table 2: Initial conditions for two stream instability and non-linear Landau damping simulations.

Param. Two-str. inst. Non-lin. Landau damp. Description

L 2π 4π Simulation domain length

T 100 60 End time

∆t 0.05 0.05 Time step

N F 199 199 Number of Fourier modes

N H 250 250 Number of Hermite modes

v th 1/ √

8 1 Thermal velocity

ν 5 1 Collisional parameter

u avg ±1 0 Average velocity

n 0.5 + 0.5 1 Average density

qom −1 −1 Charge over mass ratio

 0.001 0.5 Initial perturbation

α tol 0.1 0.01 Scaling parameter tolerance

u tol 0.01 0.01 Shift parameter tolerance

N x 500 500 # of points in real spatial space

N v 500 500 # of points in real velocity space

Two stream instability can be initiated by two oppositely-directed beams.

While the instability is growing, beam kinetic energy is transformed to the thermal energy and the potential energy of the electric field. This instability is particularly difficult to simulate with a spectral method, because the thermal and drift velocities of each species may change considerably in a time. If the Hermite basis is not optimal and not enough Hermite modes are used, then the numerical instabilities similar to those illustrated in Figures 4 and 5 may occur.

As a result the distribution function can suffer from loss of positivity or in severe cases even crash the simulation.

To initialize the simulation two electron species are used where each one

represents one beam. The initial conditions given in Table 2. Two simulations

with and without change of basis are run and results are compared.

References

Related documents

The begining of Knot theory is in a brief note made in 1833, where Carl Friedrich Gauss introduces a mathematical formula that computes the linking number of two space curves [0]..

Figure 5: These are the dierent graphs for the time of three sets of simulations for three software mod- els: one serial Fortran algorithm (square-blue-line), one threaded

To investigate if a relationship exists between price of electricity and the number of workers in the manufacturing industry, a Vector Autoregressive model test will be

Ett relativt stort antal arter registrerades dven utefter strdckor med niira an- knytning till naturbetesmarker (striickorna 5, 6.. = 9,

För att stärka sin position hade Nohab genom Gunnar Dellners försorg bildat allianser med andra industriella aktörer, något som skulle kulminera en dryg månad senare när Flygvapnets

The women would be asked to present objects that belong to them and represent for them the idea of ‘home’ and the stories related to the objects, which would explain why

We find that empirically random maps appear to model the number of periodic points of quadratic maps well, and moreover prove that the number of periodic points of random maps

We investigate the number of periodic points of certain discrete quadratic maps modulo prime numbers.. We do so by first exploring previously known results for two particular