Issues in Solving the Boltzmann Equation for Aerospace Applications - - PowerPoint PPT Presentation

Issues in Solving the Boltzmann Equation for Aerospace Applications

ICERM, Providence, RI

A high-order positivity preserving method for the Vlasov-Poisson system

James A. Rossmanith 1, David C. Seal 2, Andrew J. Christlieb 2

1Iowa State University 2Michigan State University

June 7th, 2013

Outline

Introduction Hybrid discontinous-Galerkin method Examples Other Extensions & Conclusions

Outline

Introduction Hybrid discontinous-Galerkin method Examples Other Extensions & Conclusions

Vlasov-Poisson System: A Kinetic Plasma Model

What is a plasma?

◮ ‘Fourth’ state of matter, and most abundant form of ordinary matter.

Given enough energy, anything will turn into plasma.

◮ Can think of plasma as an ionized gas or liquid. Therefore model

equations must respect E&M forces plus liquid/gas dynamics.

Scientific and Engineering Applications of Plasma

◮ Astrophysics (stars, solar corona/wind) ◮ Nuclear experiments, tocamaks, stellerators ◮ Lightning/Polar Aurorae ◮ Flourescent lights, TVs, etc.

Mathematical Models

◮ Fluid Models (e.g. continuum equations in gas dynamics) ◮ Kinetic Models - Computationally expensive.

Vlasov-Poisson System: A Collisionless Plasma

◮ Probability distribution function (PDF):

f(t, x, v) := probability of finding an electron at t & x, w/ vel. v

◮ Boltzmann equation (non-dimensional units):

∂f ∂t + v · ∇xf + a · ∇vf = C(f) Two species Vlasov-Poisson (Electrostics + C(f) ≡ 0) a = −∇φ, −∇2φ = ρ(t, x) − ρ0.

◮ Some numerical challenges

• 1. High dimensionality (3 + 3 + 1) - even worse for Boltzmann!
• 2. Conservation of mass & total energy, positivity
• 3. Complex geometries, boundary conditions.
• 4. Small time steps due to v ∈ R3 or strong electric fields.

Eulerian schemes

Selecting a scheme for Vlasov-Poisson

◮ Popular “finite-X” methods:

• 1. Finite-Difference
• 2. Finite-Volume
• 3. Finite-Element
• 4. Discontinuous-Galerkin (high-order + unstructured grids)

◮ Main advantage:

• 1. Fast convergence
• 2. Provable accuracy, load balancing for parallel computing.

◮ Main disadvantage:

• 1. Curse of dimensionality
• 2. Courrant-Fredricks-Lewy (CFL) limits

◮ M th-order discontinuous-Galerkin: ν ≈ 1/ (2M − 1).

• 3. Diffusion (in velocity space) - bad!

◮ Some recent results for Vlasov-Poisson:

[Banks & Hittinger, ’10], [Heath, Gamba, Morrison, Michler ’11], [Cheng, Gamba, Morrison ’13], . . .

Particle-in-cell and semi-Lagrangian methods

Selecting a scheme for Vlasov-Poisson

• 1. Particle in Cell
• 2. semi-Lagrangian

1.

◮ Discretize f with “macro-particles.” ◮ Key difficulties: statistical noise ∼ O(1/

√ N).

◮ [Birdsall & Langdon, 1985], [Hockney & Eastwood, 1989]

2.

◮ Main advantage: Remove CFL constraint, maintain mesh rep., . . . ◮ Main disadvantage: Boundary conditions and unstructured grids?

Mesh distortions from Lagrangian evolution + interpolation error.

◮ [Parker & Hitchon, ’97], [Sonnendrücker et al. ’99],

[Güçlü & Hitchon, ’12], . . .

Outline

Introduction Hybrid discontinous-Galerkin method Examples Other Extensions & Conclusions

Splitting Techniques for Vlasov-Poisson

Semi-Lagrangian w/ Strang splitting for Vlasov:

[Cheng & Knorr, 1976], [Besse & Sonnendrücker, ’03], [Qiu & Christlieb, ’09], . . . Split f,t + v · f,x + ∇φ(t, x) · f,v = 0 into: 1.

∆t 2 step on: f,t + v · f,x = 0

• 2. Solve
• ∇2φ = ρn+ 1

2 − ρ0

and compute En+ 1

2 = −∇φ.

• 3. ∆t step on: f,t + En+ 1

2 · f,v = 0.

4.

∆t 2 step on: f,t + v · f,x = 0

High-order Runge-Kutta Nyström splitting for V-P is an option. [Rossmanith & S, ’10], [Crouseilles et al, ’11], [S, ’12]

Multi-D extensions of SLDG?

◮ Natural extension =

⇒ Cartesian.

◮ Hybrid = Which method do we want to apply to sub-problems?

Hybrid SLDG (HSLDG) for Vlasov-Poisson

Proposed hybrid DG scheme

• 1. SLDG (Semi-Lagrangian DG) on velocity space: f,t + En+ 1

2 · f,v = 0.

[Rossmanith & S, ’10], [Einkemmer & Ostermann, ’12], not [Restelli et al, ’06], [Qiu & Shu, ’11] (less efficient).

◮ Removes CFL condition.

• 2. RKDG (Runge-Kutta DG) on configuration space: f,t + v · f,x = 0.

[Reed & Hill, 1976], [Cockburn & Shu, 90’s]

◮ Unstructured grids and Boundary conditions. ◮ Sub-cycle independent problems. ◮ Disclosure: SLDG better for structured grids!

Building blocks

Need to define how sub-problems are tackled.

Semi-Lagrangian DG (SLDG)

Advection Equation: f,t + a f,v = 0

• 1. Reconstruct. The original projection:

F (k),n

i

:= 1 ∆v

• Ci

ϕ(k)

i

(v)f(tn, v) dv; f(t, v) =

• i,k

F (k)

i

(t)ϕ(k)

i

(v).

• 2. Evolve. The evolution step is simply the exact solution:

f(t, v) = f0(v − at); f0(v) = f(0, v).

• 3. Average. This step requires the solution from step 2.

F (k),n+1

i

:= 1 ∆v

• Ci

ϕ(k)

i

(v)f(tn+1, v) dv = 1 ∆v

• Ci

ϕ(k)

i

(v)f(tn, v − a(t − tn))) dv

2D Semi-Lagrangian DG

Constant coefficient advection: a.k.a. corner transport + high-order

• 1. Forward (cell discontinuities)
• 2. Backward (quadrature points)

Advection equation: f,t − E1 f,vx − E2f,vy = 0

F (k),n+1

i,j

:= 1

• Ti,j
• Ti,j

ϕ(k)

i,j (x, y) f h

tn+1, x, v

• dx dv

= 1

• Ti,j
• 4
• m=1
• T m

i,j

ϕ(k)

i,j (x, y)f h

tn, x + ∆tE1, v + ∆tE2 dx dv

◮ Mass conservation (and stability!) come from exact integration

Hybrid DG: Reduction to sub-problems

Not trivial because basis functions depend on transverse variable

• 1. Reconstruct: Start with a full 2N-dimensional solution:

f h (tn, x, v)

• T 2N

i

=

• k2

F (k2)

i

(tn)ϕ(k2)

2N (x, v)

• 2. Project full solution onto to sub-problems (N-dimensional) by taking

slices at quadrature points {v1, v2, . . . vMN }: f h (tn, x, vm)

• T N

i

=

• k

F (k)

N,i(tn)ϕ(k) N (x)

• 3. Evolve (RKDG or SLDG) sub-problems:

f h tn+1, x, vm

• T N

i

=

• k

F (k)

N,i(tn+1)ϕ(k) N (x)

• 4. Integrate: Sub-problems integrated up to 2N-dimensional problem:

f h tn+1, x, v

• T 2N

i

=

• k2

F (k2)

i

(tn+1)ϕ(k2)

2N (x, v)

1D-1V Example: f,t + v f,x = 0

• 1. Full 2D-solution.
• 2. 1D-Problems (at quadrature points).
• 3. Evolve lower-D problems (SLDC/RKDG)
• 4. Integrate up to full soln.

HSLDG for Vlasov-Poisson: Sub-Cycling

Basic Idea: March sub-problems at their own pace

◮ Many independent sub-Problems: f,t + vmf,x = 0. Global restriction:

∆t ≤ CFL ∆x |vmax| ≪ CFL ∆x |vm| ∼ ∆tlocal

◮ For each equation, define: ∆tlocal := ∆t/ max

• 1,
• |vm|∆t

ν∆x

• .

Outline

Introduction Hybrid discontinous-Galerkin method Examples Other Extensions & Conclusions

Example: Forced Vlasov-Poisson Test Problem

◮ Method of Manufactured Solutions produces a source term:

f,t + vf,x + E(t, x)f,v = ψ(t, x, v), E,x = ∞

−∞

f(t, x, v) dv − √π.

◮ Choose the exact solution: f(t, x, v) = {2 − cos (2x − 2πt)} e− 1

4 (4v−1)2.

◮ New Hybrid Operators:

Sub-cycled RKDG on Problem A: f,t + v f,x = ψ(t, x, v), SLDG on Problem B: f,t + E(t, x) f,v = 0,

Mesh HSL2 Error log2(Ratio) HSL4 Error log2(Ratio)

102 5.193 × 10−1 – 1.085 × 100 – 202 1.432 × 10−1 1.858 2.725 × 10−1 1.993 402 1.640 × 10−2 3.126 1.652 × 10−2 4.044 802 3.438 × 10−3 2.255 7.058 × 10−4 4.548 1602 8.333 × 10−4 2.045 3.421 × 10−5 4.367 3202 2.068 × 10−4 2.011 1.953 × 10−6 4.131 6402 5.161 × 10−5 2.003 1.197 × 10−7 4.028

Example: Landau Damping with HSLDG

Initial Conditions: f(t = 0, x, v) =

• 1 + α cos(kx)
• 1

√ 2π e− v2

2

Weak Landau Strong Landau

10 20 30 40 50 60 10

−8

10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10 20 30 40 50 60 10

−4

10

−3

10

−2

10

−1

10 10

1

α = 0.01 α = 0.5

Plasma Sheath Simulation (1D-1V) with HSLDG

Initial conditions (in dimensional quantitites): ˜ f(0, ˜ x, ˜ v) = ˜ f0(˜ v) = ˜ ρ0

• 2πmk˜

θ0 exp

• − m˜

v2 2k˜ θ0

• ,

˜ x ∈ [0, L]. Boundary conditions: ˜ f(t, ˜ x = 0, ˜ v) = 0, ˜ f(t, ˜ x = L, ˜ v) = ˜ f0(˜ v), ˜ φ(t, ˜ x = 0) = 0, ˜ φ˜

x(t, ˜

x = L) = 0. 1 eV electrons, 0.1m domain, density ˜ ρ0 = 9.10938188 × 10−18 kg

m .

Plasma Sheath Simulation (1D-1V) with HSLDG

(sheath-movie.mp4)

Plasma Sheath Simulation (1D-1V) with HSLDG

Plasma Sheath Simulation (1D-1V) with HSLDG

0.02 0.04 0.06 0.08 0.1 2.5 5 7.5 10 12.5 x 10

12

ne(t,x) at t = 3.0269e−06 0.02 0.04 0.06 0.08 0.1 −800 −600 −400 −200 200 400 600 800 E(t,x) at t = 3.0269e−06

Number density Electric field

Plasma Sheath Simulation (2D-2V) with HSLDG

◮ Initial conditions: constant density, ˜

ρ0 = 9.10938188 × 10−18 kg

m2 .

◮ Boundary conditions: ˜

f(˜ t, ˜ x = L, ˜ v) = 0, ˜ φ(˜ t, ˜ x = L) = 0.

−0.1 −0.05 0.05 0.1 −0.1 −0.08 −0.06 −0.04 −0.02 0.02 0.04 0.06 0.08 0.1 rho(t,x,y) at t = 0 [DoGPack]

0.02 0.04 0.06 0.08 0.1 2.5 5 7.5 10 12.5 x 10

12

ne(t,r) at t = 0 [DoGPack]

1 eV electrons,

• π · 0.12

m2 domain.

Plasma Sheath Simulation (2D-2V) with HSLDG

−0.1 −0.05 0.05 0.1 −0.1 −0.08 −0.06 −0.04 −0.02 0.02 0.04 0.06 0.08 0.1 rho(t,x,y) at t = 4.4843e−08 [DoGPack]

0.02 0.04 0.06 0.08 0.1 2.5 5 7.5 10 12.5 x 10

12

ne(t,r) at t = 4.4843e−08 [DoGPack]

Plasma Sheath Simulation (2D-2V) with HSLDG

−0.1 −0.05 0.05 0.1 −0.1 −0.08 −0.06 −0.04 −0.02 0.02 0.04 0.06 0.08 0.1 rho(t,x,y) at t = 6.7265e−08 [DoGPack]

0.02 0.04 0.06 0.08 0.1 2.5 5 7.5 10 12.5 x 10

12

ne(t,r) at t = 6.7265e−08 [DoGPack]

Plasma Sheath Simulation (2D-2V) with HSLDG

−0.1 −0.05 0.05 0.1 −0.1 −0.08 −0.06 −0.04 −0.02 0.02 0.04 0.06 0.08 0.1 rho(t,x,y) at t = 1.1211e−07 [DoGPack]

0.02 0.04 0.06 0.08 0.1 2.5 5 7.5 10 12.5 x 10

12

ne(t,r) at t = 1.1211e−07 [DoGPack]

Plasma Sheath Simulation (2D-2V) with HSLDG

−0.1 −0.05 0.05 0.1 −0.1 −0.08 −0.06 −0.04 −0.02 0.02 0.04 0.06 0.08 0.1 rho(t,x,y) at t = 2.2422e−07 [DoGPack]

0.02 0.04 0.06 0.08 0.1 2.5 5 7.5 10 12.5 x 10

12

ne(t,r) at t = 2.2422e−07 [DoGPack]

Plasma Sheath Simulation (2D-2V) with HSLDG

−0.1 −0.05 0.05 0.1 −0.1 −0.08 −0.06 −0.04 −0.02 0.02 0.04 0.06 0.08 0.1 rho(t,x,y) at t = 5.6054e−07 [DoGPack]

0.02 0.04 0.06 0.08 0.1 2.5 5 7.5 10 12.5 x 10

12

ne(t,r) at t = 5.6054e−07 [DoGPack]

Plasma Sheath Simulation (2D-2V) with HSLDG

−0.1 −0.05 0.05 0.1 −0.1 −0.08 −0.06 −0.04 −0.02 0.02 0.04 0.06 0.08 0.1 rho(t,x,y) at t = 6.7265e−07 [DoGPack]

0.02 0.04 0.06 0.08 0.1 2.5 5 7.5 10 12.5 x 10

12

ne(t,r) at t = 6.7265e−07 [DoGPack]

Plasma Sheath Simulation (2D-2V) with HSLDG

Cross sections of velocity space

−0.2 −0.1 0.1 0.2 −0.2 −0.15 −0.1 −0.05 0.05 0.1 0.15 0.2 f(t,vx,vy) at t = 0 [DoGPack]

50 100 150 200 250 300

−0.2 −0.1 0.1 0.2 −0.2 −0.15 −0.1 −0.05 0.05 0.1 0.15 0.2 f(t,vx,vy) at t = 0 [DoGPack]

−25 −20 −15 −10 −5 5 10 15 20 25

Plasma Sheath Simulation (2D-2V) with HSLDG

Cross section of velocity space

−0.2 −0.1 0.1 0.2 −0.2 −0.15 −0.1 −0.05 0.05 0.1 0.15 0.2 f(t,vx,vy) at t = 0 [DoGPack]

−25 −20 −15 −10 −5 5 10 15 20 25

−0.2 −0.1 0.1 0.2 −0.2 −0.15 −0.1 −0.05 0.05 0.1 0.15 0.2 f(t,vx,vy) at t = 5.6054e−08 [DoGPack]

−25 −20 −15 −10 −5 5 10 15 20 25

−0.2 −0.1 0.1 0.2 −0.2 −0.15 −0.1 −0.05 0.05 0.1 0.15 0.2 f(t,vx,vy) at t = 5.6054e−07 [DoGPack]

−25 −20 −15 −10 −5 5 10 15 20 25

−0.2 −0.1 0.1 0.2 −0.2 −0.15 −0.1 −0.05 0.05 0.1 0.15 0.2 f(t,vx,vy) at t = 6.7265e−07 [DoGPack]

−25 −20 −15 −10 −5 5 10 15 20 25

Outline

Introduction Hybrid discontinous-Galerkin method Examples Other Extensions & Conclusions

Conclusions & Future Work

◮ Vlasov-Poisson: nonlinear/nonlocal advection equation ◮ Semi-Lagrangian and Eulerian grid-based methods as alternative

to PIC and particle methods

◮ Operator splitting + sub-cycling + semi-Lagrangian relaxes CFL. ◮ 1D-1V and preliminary 2D-2V results: collisionless simulation of

the formation of a plasma sheath, Landau damping, bump on tail, ... http://www.dogpack-code.org

◮ The Hybrid HSLDG and the SLDG Method attain:

• 1. Unconditionally stable; (sub-cycling / semi-Lagrangian)
• 2. High-order space (5th order) and time (4th order);
• 3. Mass conservative; and
• 4. Positivity-preserving

Future work

◮ Extra options for 2D-Poisson solver. ◮ Vlasov-Maxwell, SLDG on even higher dimensions, . . .

Conclusions Thank You!

This work was partially supported by the following agencies:

◮ Michigan State University Foundation SPG-RG100059 ◮ Air Force Office of Scientific Research (AFOSR)

FA9550-11-1-0281, FA9550-12-1-0343, FA9550-12-1-0455

◮ National Science Foundation (NSF) DMS-1115709

The DoGPack software package

http://www.dogpack-code.org

Semi-Lagrangian discontinuous Galerkin

Positivity Preserving Limiter: Step I Provided cell average is positive, one can gaurantee as many points as requested are positive without losing order of accuracy.

• 1. To limit cell Tij, compute

m := min

(x,v)∈Tij

f h(x, v),

• r at least an approximate minimum.
• 2. Define the limited solution [Zhang & Shu, 2010] as

˜ f h(ξ, η) := F (1) + θ ·

• f h(ξ, η) − F (1)

, θ = min

• 1,

F (1) F (1) − m

• ,

0 ≤ θ ≤ 1.

• 3. Main advantages: Locally applied limiter, easy to implement.
• 4. Difficulties:

◮ Exact minimum of the function on each cell? ◮ How does one make sure average is positive?

Semi-Lagrangian discontinuous Galerkin

Positivity Preserving Limiter [Rossmanith and S., 2011]: Step II Theorem: A single advect-project step preserves positivity in the mean, and hence the

• verall algorithm is positivity preserving. That is,
• F (1)

i,j

n+1 := 1

• Ti,j
• Ti,j

f h tn+1, x, v

• dx dv ≥ 0.

Proof:

◮ The points chosen for integration are forced to be positive via a

judicious and finite choice of limiting points in a preprocessing step:

◮ Step I makes solution positive at our choice of quadrature points. ◮ All quadrature weights are positive. ◮ Extra effort reduces number of points to 2MK, where K :=

M

2

• .

Semi-Lagrangian DG-FEM

4th order-in-time

◮ Poisson’s equation

−∇2φ = ρn − ρ0

◮ Time derivatives

En = ∇φ, En

,t = − (ρu)n ,

En

,tt = ∇ · En − (ρE)n

En

,ttt = ∇ · (ρuE + ρEu)n − ∇ · ∇ · Fn + En ∇ · (ρu)n +

• ρ2u

n ¯ E(t, x) = En + (t − tn) En

,t + 1

2(t − tn)2 En

,tt + 1

6(t − tn)3 En

,ttt

Fourth-order splitting:

Stage 1: + 0.6756 ∆t step on f,t + v · f,x = 0 Stage 2: + 1.3512 ∆t step on f,t + ¯ E (t, x) · f,v = 0 Stage 3: − 0.1756 ∆t step on f,t + v · f,x = 0 Stage 4: − 1.7024 ∆t step on f,t + ¯ E (t, x) · f,v = 0 Stage 5: − 0.1756 ∆t step on f,t + v · f,x = 0 Stage 6: + 1.3512 ∆t step on f,t + ¯ E (t, x) · f,v = 0 Stage 7: + 0.6756 ∆t step on f,t + v · f,x = 0

Eulerian Discontinuous Galerkin methods

Spatial discretization [Cockburn & Shu, 1990’s]

◮ Basis functions:

ϕ(ℓ)(x) =

• 1, x, y, x2, xy, y2, . . .
• ◮ Galerkin expansion:

f h(t, x) = M(M+1)/2

k=1

F (ℓ)(t) ϕ(ℓ)(x)

◮ Semi-discrete weak-form of f,t + ∇ · R(f) = 0:

• T

ϕ(ℓ) f,t dx = −

• T

ϕ(ℓ) ∇ · R(f) dx = ⇒ d dtF (ℓ) = 1 |T |

• T

∇ϕ(ℓ) · R(f) dx

• Interior

− 1 |T |

• ∂T

ϕ(ℓ) R(f) · ds

• Edge

◮ Interior: quadrature,

Edge: quadrature, then approx Riemann soln