Implicit Methods and the M3D- C 1 Approach Stephen C. Jardin - - PowerPoint PPT Presentation

MHD Simulations for Fusion Applications

Lecture 4

Implicit Methods and the M3D-C1 Approach

Stephen C. Jardin Princeton Plasma Physics Laboratory

1

CEMRACS ‘10 Marseille, France July 22, 2010

Outline of Remainder of Lecture

• Galerkin method
• C1 finite elements
• Examples
• Split implicit time

differencing

Yesterday Today

• Split implicit time

differencing for MHD

• Accuracy, spectral

pollution, and representation of the vector fields

• Projections of the split

implicit equations

• Energy conserving subsets
• 3D solver strategy
• Examples

3

2-Fluid MHD Equations:

( )

2 2

1 ( ) continuity Maxwell ( ) momentum Ohm's law 3 3 electron energy 2 2 3 3 2 2

e i e e e G i e F i V e

p ne n n t t nM p t p p p t p J p t Q S μ η η μ

Δ

∂ + ∇ • = ∂ ∂ = −∇× ∇ = = ∇× ∂ ∂ +

• ∇

+ ∇ = × − ∂ + × = + ∂ ⎛ ⎞ + ∇ = − ∇ × − ∇ ⎜ ⎟ ∂ + ⎝ ⎠ ∂ ⎛ ⎞ + ∇ = ⎜ ⎟ ∂ ⎝ ⎠ ∇ ∇ + −∇ + + V B E B J B V V V J B E V B Π J B V V J q V V i i i i i i

2

ion energy

i Fi i

S p V Q μ

Δ

+ ∇ −∇ − + − ∇ q V i i

Resistive MH 2- Idea flui D d l MHD MHD

number density magnetic field current density electric field mass density

i

n nM ρ ≡ Β J E fluid velocity electron pressure ion pressure electron charge

e i e i

p p p p p e ≡ + V viscosity resistivity heat fluxes equipartition permeability Q μ μ η

Δ i e

q ,q

4

The split-implicit time advance

5

1 1 1 1/2 1/2 1/2 1/2 1 1 1 1/2 1/2 1 1

(1 ) (1 )

n n n n n n j j j j j j n n n n n n j j j j j j

u u v v v v c t x x v v u u u u c t x x θ θ δ δ δ θ θ δ δ δ

+ + + + − + − + + + + + + +

⎡ ⎤ ⎛ ⎞ ⎛ ⎞ − − − = + − ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎝ ⎠ ⎣ ⎦ ⎡ ⎤ ⎛ ⎞ ⎛ ⎞ − − − = + − ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎝ ⎠ ⎣ ⎦

( ) ( )

1 1 1 1 1 1 1 1/2 1/2 1 2 2 2 2 1 1 1 1/2 1/2 1 1

2 2 ( ) (1 ) (1 )

n n n n n n n n j j j j j j j j n n j j n n n n n n j j j j j j

u u u u u u v v u u tc tc x x x tc v v u u u u x δ θ θ θ δ δ δ δ δ θ θ δ

+ + + + − + − + − + + + + + + + +

⎡ ⎤ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ − + − + − = + + − + ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎣ ⎦ ⎡ ⎤ = + − + − − ⎣ ⎦

Consider a simple 1-D Hyperbolic System of Equations (Wave Equation) Substitute from second equation into first:

These two equations are completely equivalent to those on the previous page, but can be solved sequentially! Only first involves Matrix Inversion … Diagonally Dominant

6

u v c t x v u c t x ∂ ∂ = ∂ ∂ ∂ ∂ = ∂ ∂

n n

u v c v t t x t v u c u t t x t θδ θδ ∂ ∂ ∂ ⎡ ⎤ = + ⎢ ⎥ ∂ ∂ ∂ ⎣ ⎦ ∂ ∂ ∂ ⎡ ⎤ = + ⎢ ⎥ ∂ ∂ ∂ ⎣ ⎦

An alternate derivation:

Expand RHS in Taylor series in time to time- center

7

u v c t x v u c t x ∂ ∂ = ∂ ∂ ∂ ∂ = ∂ ∂

n n

u v c v t t x t v u c u t t x t θδ θδ ∂ ∂ ∂ ⎡ ⎤ = + ⎢ ⎥ ∂ ∂ ∂ ⎣ ⎦ ∂ ∂ ∂ ⎡ ⎤ = + ⎢ ⎥ ∂ ∂ ∂ ⎣ ⎦

n n n

u u c v t c u t t x x t v u c u t t x t θδ θδ θδ ⎡ ⎤ ∂ ∂ ⎛ ∂ ∂ ⎞ ⎡ ⎤ = + + ⎢ ⎥ ⎜ ⎟ ⎢ ⎥ ∂ ∂ ∂ ∂ ⎣ ⎦ ⎝ ⎠ ⎣ ⎦ ∂ ∂ ∂ ⎡ ⎤ = + ⎢ ⎥ ∂ ∂ ∂ ⎣ ⎦

An alternate derivation:

Expand RHS in Taylor series in time to time- center Substitute from second equation into first

8

2 2 2 2 2 1 2 2 2 2 1 1

1 ( ) 1 (1 )( ) (1 )

n n n n n n n

t c u t c u tc v x x x v v tc u u x x θ δ θ θ δ δ δ θ θ

+ + +

⎡ ⎤ ⎡ ⎤ ∂ ∂ ∂ − = + − + ⎢ ⎥ ⎢ ⎥ ∂ ∂ ∂ ⎣ ⎦ ⎣ ⎦ ∂ ∂ ⎡ ⎤ = + + − ⎢ ⎥ ∂ ∂ ⎣ ⎦ u v c t x v u c t x ∂ ∂ = ∂ ∂ ∂ ∂ = ∂ ∂

n n

u v c v t t x t v u c u t t x t θδ θδ ∂ ∂ ∂ ⎡ ⎤ = + ⎢ ⎥ ∂ ∂ ∂ ⎣ ⎦ ∂ ∂ ∂ ⎡ ⎤ = + ⎢ ⎥ ∂ ∂ ∂ ⎣ ⎦

n n n

u u c v t c u t t x x t v u c u t t x t θδ θδ θδ ⎡ ⎤ ∂ ∂ ⎛ ∂ ∂ ⎞ ⎡ ⎤ = + + ⎢ ⎥ ⎜ ⎟ ⎢ ⎥ ∂ ∂ ∂ ∂ ⎣ ⎦ ⎝ ⎠ ⎣ ⎦ ∂ ∂ ∂ ⎡ ⎤ = + ⎢ ⎥ ∂ ∂ ∂ ⎣ ⎦

An alternate derivation:

Expand RHS in Taylor series in time to time- center Substitute from second equation into first

Use standard centered difference in time:

9

u v c t x v u c t x ∂ ∂ = ∂ ∂ ∂ ∂ = ∂ ∂

n n

u v c v t t x t v u c u t t x t θδ θδ ∂ ∂ ∂ ⎡ ⎤ = + ⎢ ⎥ ∂ ∂ ∂ ⎣ ⎦ ∂ ∂ ∂ ⎡ ⎤ = + ⎢ ⎥ ∂ ∂ ∂ ⎣ ⎦

n n n

u u c v t c u t t x x t v u c u t t x t θδ θδ θδ ⎡ ⎤ ∂ ∂ ⎛ ∂ ∂ ⎞ ⎡ ⎤ = + + ⎢ ⎥ ⎜ ⎟ ⎢ ⎥ ∂ ∂ ∂ ∂ ⎣ ⎦ ⎝ ⎠ ⎣ ⎦ ∂ ∂ ∂ ⎡ ⎤ = + ⎢ ⎥ ∂ ∂ ∂ ⎣ ⎦

2 2 2 2 2 1 2 2 2 2 1 1

1 ( ) 1 (1 )( ) (1 )

n n n n n n n

t c u t c u tc v x x x v v tc u u x x θ δ θ θ δ δ δ θ θ

+ + +

⎡ ⎤ ⎡ ⎤ ∂ ∂ ∂ − = + − + ⎢ ⎥ ⎢ ⎥ ∂ ∂ ∂ ⎣ ⎦ ⎣ ⎦ ∂ ∂ ⎡ ⎤ = + + − ⎢ ⎥ ∂ ∂ ⎣ ⎦

An alternate derivation:

Expand RHS in Taylor series in time to time- center Substitute from second equation into first

Use standard centered difference in time:

This is the same operator as before when centered spatial differences used

W δ

(next 3 vgs)

let = introduce a displacement vector linearize equations about =0 t ∂ ∂ ξ V ξ V

( ) ( ) ( ) ( )

2 1 1 1 2 1 1

1 1 2 3 p t p p p ρ μ μ ∂ + ∇ = ∇× × + ∇× × ∂ = ∇× × + ∇ = − ∇ ξ B B B B B ξ B ξ ξ i i

( ) ( )

( )

( )

5 3 2 2

• r,

Linearized equation of motion 1 ( ) ( ) t p p μ ρ ⎡ ⎤ = ∇× × + ∇× × + ∇ ∇ + ∇ ⎣ ∂ = ⎦ ≡ = ∇× × ∂

1

ξ F F ξ B Q Q B ξ ξ Q B ξ ξ B i i

Linear ideal MHD response can be analyzed with δW approach

Ideal MHD

• perator

( ) ( )

( )

( )

2 2 5 3

1 ( ) ( ) p p t ρ μ ⎡ ⎤ = ∇× × + ∇× × + ∇ ∇ + ∇ ⎣ ⎦ ≡ = ∇ = ∂ × × ∂

1

ξ F ξ B Q Q B ξ ξ Q B B F ξ ξ i i

† 1 2

Assume ( , ) ( ) and take dot product with - , and integrate over volume to get perturbed energy

i t

x t x e ω = ξ ξ ξ

2 2 † 1 1 2 2 †

( , ) if 0 for any displacement field instabilit ( y ) d d W W ρω τ τ δ δ = − < → ≡

∫ ∫

ξ ξ ξ ξ ξ F ξ i

δW is the potential energy for a given displacement field

[ ]

( )( )

2 2 5 1 3 2 1 1 2 2

2 2 p d p B W

μ μ

δ τ σ

⊥ ⊥ ⊥ ⊥ ⊥ ⊥ ⊥

∇ ⎡ ⎤ + + ∇ ⎢ ⎥ = ⎢ ⎥ − ∇ − × ⎣ ⎦ +

∫

Q ξ ξ κ ξ B Q κ ξ ξ ξ i i i i i i

( )

2

/ unit vector in direction of field curvature of magnetic field / parallel current density perturbed magnetic field perpendicular component of perpendicular component of B B σ

⊥ ⊥

= = ∇ = = ∇× × = = b B κ b b J B Q ξ B Q Q ξ ξ i i

The plasma motion will be such as to minimize δW

This is the term associated with the fast wave…it is positive definite with a large multiplier. Any unstable plasma motion will make this term small

2 ~ 0

⊥ ⊥

∇ ≅ − ξ ξ κ i i

14

[ ] [ ]

1 p p p p ρ μ γ = ∇× × −∇ = ∇× × = − ∇ − ∇ V B B B V B V V

• i

i

Now apply the split implicit advance to the basic 3D (ideal) MHD equations:

Ideal MHD Equations for velocity, magnetic field, and pressure: Symmetric Hyperbolic System 7-waves

15

[ ] [ ]

1 p p p p ρ μ γ = ∇× × −∇ = ∇× × = − ∇ − ∇ V B B B V B V V

• i

i

( ) ( )

( )

( ) ( ) ( )

1 t t p tp t p t p p t ρ θδ θδ θδ μ θδ θδ γ θδ ⎡ ⎤ = ∇× + × + −∇ + ⎣ ⎦ ⎡ ⎤ = ∇× + × ⎣ ⎦ = − + ∇ − ∇ + V B B B B B V V B V V V V

• i

i

Ideal MHD Equations for velocity, magnetic field, and pressure: Symmetric Hyperbolic System 7-waves Taylor Expand in Time as before

Now apply the split implicit advance to the basic 3D (ideal) MHD equations:

16

[ ] [ ]

1 p p p p ρ μ γ = ∇× × −∇ = ∇× × = − ∇ − ∇ V B B B V B V V

• i

i

( ) ( )

( )

( ) ( ) ( )

1 t t p tp t p t p p t ρ θδ θδ θδ μ θδ θδ γ θδ ⎡ ⎤ = ∇× + × + −∇ + ⎣ ⎦ ⎡ ⎤ = ∇× + × ⎣ ⎦ = − + ∇ − ∇ + V B B B B B V V B V V V V

• i

i

{ } { }

2 2 1 2

1 ( ) ( 1)( ) ( )

n n n

t L t L t p ρ θ δ ρ θ θ δ δ μ

+

⎧ ⎫ − = − − + −∇ + ∇× × ⎨ ⎬ ⎩ ⎭ V V B B

{ } ( )

{ }

( ) ( ) ( )

1 1 L p p μ μ γ = ∇× ∇× × × + ∇× × ∇× × ⎡ ⎤ ⎡ ⎤ ⎣ ⎦ ⎣ ⎦ +∇ ∇ + ∇ V V B B B V B V V i i

Ideal MHD Equations for velocity, magnetic field, and pressure: Symmetric Hyperbolic System 7-waves Taylor Expand in Time as before Substitute from 2nd and 3rd equation into first, finite difference in time: MHD Operator:

Now apply the split implicit advance to the basic 3D (ideal) MHD equations:

17

[ ] [ ]

1 p p p p ρ μ γ = ∇× × −∇ = ∇× × = − ∇ − ∇ V B B B V B V V

• i

i

( ) ( )

( )

( ) ( ) ( )

1 t t p tp t p t p p t ρ θδ θδ θδ μ θδ θδ γ θδ ⎡ ⎤ = ∇× + × + −∇ + ⎣ ⎦ ⎡ ⎤ = ∇× + × ⎣ ⎦ = − + ∇ − ∇ + V B B B B B V V B V V V V

• i

i

{ } { }

1/2 2 2 1 2 2

1 ( ) ( ) ( )

n n n

t L t L t p ρ θ δ ρ θ δ δ μ

+ +

⎧ ⎫ − = − + −∇ + ∇× × ⎨ ⎬ ⎩ ⎭ V V B B

{ } ( )

{ }

( ) ( ) ( )

1 1 L p p μ μ γ = ∇× ∇× × × + ∇× × ∇× × ⎡ ⎤ ⎡ ⎤ ⎣ ⎦ ⎣ ⎦ +∇ ∇ + ∇ V V B B B V B V V i i

Now apply this technique to the basic 3D MHD equations: Ideal MHD Equations for velocity, magnetic field, and pressure: Symmetric Hyperbolic System 7-waves Taylor Expand in Time as before Substitute from 2nd and 3rd equation into first, finite difference in time: MHD Operator: note!

18

{ } { }

1/2 2 2 1 2 2

1 ( ) ( ) ( )

n n n

t L t L t p ρ θ δ ρ θ δ δ μ

+ +

⎧ ⎫ − = − + −∇ + ∇× × ⎨ ⎬ ⎩ ⎭ V V B B

{ } ( )

{ }

( ) ( ) ( )

1 1 L p p μ μ γ = ∇× ∇× × × + ∇× × ∇× × ⎡ ⎤ ⎡ ⎤ ⎣ ⎦ ⎣ ⎦ +∇ ∇ + ∇ V V B B B V B V V i i

MHD Operator: note!

{ } { }

2 2 1 2

1 ( ) ( 1)( ) ( )

n n n

t L t L t p ρ θ δ ρ θ θ δ δ μ

+

⎧ ⎫ − = − − + −∇ + ∇× × ⎨ ⎬ ⎩ ⎭ V V B B

Advantages of

• ver:
• 1. Gives correct steady-state physics (when Vn+1 = Vn)
• 2. Gives stable numerical method for θ ≥ 1/2 even when

plasma is unstable (L has negative eigenvalues)

• 3. Still second order accurate in time (obtained by

evaluating p and B at the half-time level)

• 4. Effectively introduces a k-dependent mass term to

provide numerical stability

19

{ } { }

1/2 2 2 1 2 2

1 ( ) ( ) ( )

n n n

t L t L t p ρ θ δ ρ θ δ δ μ

+ +

⎧ ⎫ − = − + −∇ + ∇× × ⎨ ⎬ ⎩ ⎭ V V B B

{ } ( )

{ }

( ) ( ) ( )

1 1 L p p μ μ γ = ∇× ∇× × × + ∇× × ∇× × ⎡ ⎤ ⎡ ⎤ ⎣ ⎦ ⎣ ⎦ +∇ ∇ + ∇ V V B B B V B V V i i

Summary of Split Implicit Time Advance

( ) ( ) ( ) ( ) ( ) ( ) ( )

3/2 1/2 1 3/2 1/2 3/2 1/2 1 1 3/2 1/2 3/2 1/2 1/2 1/2 3/2 1/2 3/2 1/2 1/2 3/

(1 ) (1 ) (1 ) (1 ) (1 2 ) 1

n n n n n n n n n n n n n n n n n n n n n

p p t p p p p t tS t ne δ θ θ γ θ θ δ δ θ θ η θ θ θ δ θ θ

+ + + + + + + + + + + + + + + + + + + + +

= − ∇ + − ∇ − + − ∇ + × + − − ∇× + − ∇× − ∇× × = + ∇× − + ∇× × + ∇× × V V V B B B B B B B B B B B B i i

( )

2 3/2 1/2

1 (1 )

n n e e

p p ne θ θ

+ +

⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎢ ⎥ ⎢ ⎥ − ∇ + − ⎢ ⎥ ⎣ ⎦

1

• nly

implic ! it in

n+

V

3/2

• nly

implicit i ! n

n

p +

3/2

• nly

implici ! t in

n+

B

Hall terms

20

{ } { }

1/2 2 2 1 2 2

1 ( ) ( ) ( )

n n n

t L t L t p ρ θ δ ρ θ δ δ μ

+ +

⎧ ⎫ − = − + −∇ + ∇× × ⎨ ⎬ ⎩ ⎭ V V B B

{ } ( )

{ }

( ) ( ) ( )

1 1 L p p μ μ γ = ∇× ∇× × × + ∇× × ∇× × ⎡ ⎤ ⎡ ⎤ ⎣ ⎦ ⎣ ⎦ +∇ ∇ + ∇ V V B B B V B V V i i

ϕ Z R

• Need to solve this in 3D torus with

strong magnetic field in toroidal direction (ϕ) … anisotropy

• Wide range of wave speeds with

differing polarizations leads to ill- conditioned matrices

• Gradients in (R,Z) plane much larger

than in ϕ direction

• Also need to preserve

∇ = B i

The Major Challenge in M3D-C1 is in solving the Implicit Velocity Equation

21

Accuracy and Spectral Pollution

Because the externally imposed toroidal field in a tokamak is very strong, any plasma instability will slip through this field and not compress it. We need to be able to model this motion very accurately because of the weak forces causing the instability.

2 2 2 2 2

1 ( ) U R R R F f R f φ φ φ φ φ ω χ ψ

⊥

= ∇ ×∇ + ∇ + ∇ ∂ = ∇ ×∇ −∇ + + ∇ ∇ ∂ V B

( )

( )

[ ]

2

t R U F F U φ φ φ φ φ ∂ ⎡ ⎤ ∇ = ∇× × ⎢ ⎥ ∂ ⎣ ⎦ ⎡ ⎤ = ∇ ∇× ∇ ×∇ × ∇ ⎣ ⎦ = − ∇ ∇ ×∇ = B V B i i i

In M3D-C1, we express the velocity and magnetic fields as shown: Consider now the action of the first term in V on the external toroidal field: The unstable mode will mostly consist of the velocity component U. The velocity field U does not compress the external toroidal field!

big! 2

F R U φ φ = ∇ = ∇ ×∇ B V

[ ]

( )( )

2 2 5 1 3 2 1 1 2 2

2 2 p d p B W

μ μ

δ τ σ

⊥ ⊥ ⊥ ⊥ ⊥ ⊥ ⊥

∇ ⎡ ⎤ + + ∇ ⎢ ⎥ = ⎢ ⎥ − ∇ − × ⎣ ⎦ +

∫

Q ξ ξ κ ξ B Q κ ξ ξ ξ i i i i i i

( )

2

/ unit vector in direction of field curvature of magnetic field / parallel current density perturbed magnetic field perpendicular component of perpendicular component of B B σ

⊥ ⊥

= = ∇ = = ∇× × = = b B κ b b J B Q ξ B Q Q ξ ξ i i

The plasma motion will be such as to minimize δW

This is the term associated with the fast wave…it is positive definite with a large multiplier. Any unstable plasma motion will make this term small

2 ~ 0

⊥ ⊥

∇ ≅ − ξ ξ κ i i

23

{ } ( )

{ }

( ) ( ) ( )

1 1 L p p γ μ μ = ∇× ∇× × × + ∇× × ∇× × +∇ ∇ + ∇ ⎡ ⎤ ⎡ ⎤ ⎣ ⎦ ⎣ ⎦ V V B B B V B V V i i

ϕ Z R

2 2 2

1 R R U R ϕ ϕ ω χ

⊥

= ∇ ×∇ + ∇ + ∇ V

2 2 2

1 R R R U ω ϕ ϕ χ

⊥

= ∇ ×∇ + ∇ + ∇ V

• 2

11 12 13 21 22 23 31 32 33

δW( , ) ( ) δW ( ) δW ( , ) δW ( , ) δW ( , ) δW ( ) δW ( , ) δW ( , ) δW ( , ) W ) , ( , δ , d U U U R U U L U χ χ ω ω ω ω ω ω χ χ χ χ ≡ = + + + + + + + +

∫∫

V V V V

• i

This is the ideal MHD operator of Bernstein, Freeman, Kruskal, and Kulsrud (1958) Define now 2 displacement (velocity) fields: consider the functional: can be broken up into these 9 parts, each of which is a quadratic functional

24

{ } { }

1/2 2 2 1 2 2

1 ( ) ( ) ( )

n n n

t L t L t p ρ θ δ ρ θ δ δ μ

+ +

⎧ ⎫ − = − + −∇ + ∇× × ⎨ ⎬ ⎩ ⎭ V V B B

ϕ Z R

• To solve this by the finite element

method, we need to take projections to get scalar equations, and then to take the weak form of those equations.

2 2 2 2 2 2

1

i i i

d R R d R R d R R ν ϕ ν ϕ ν

⊥ ⊥

− ∇ ∇ × ∇ − ∇

∫∫ ∫∫ ∫∫

i i i

2 2 2 2 2 2

1

i i i

d R R d R R d R R ν ϕ ν ϕ ν

⊥ ⊥

∇ ×∇ ∇ ∇

∫∫ ∫∫ ∫∫

i i i

by parts

2 2 2

1 R R U R ϕ ϕ ω χ

⊥

= ∇ ×∇ + ∇ + ∇ V

( , ) is i finite element trial function

th i R Z

ν

Projection or annihilation operators: Same form as velocity!

25

{ } { }

1/2 2 2 1 2 2

1 ( ) ( ) ( )

n n n

t L t L t p ρ θ δ ρ θ δ δ μ

+ +

⎧ ⎫ − = − + −∇ + ∇× × ⎨ ⎬ ⎩ ⎭ V V B B

ϕ Z R

• Consider the effect of these projection
• perators on the MHD operator

{ } { } { }

2 2 2 2 2 2

1

i i i

d R R L d R R L d R L R ν ϕ ν ϕ ν

⊥ ⊥

∇ ×∇ ∇ ∇

∫∫ ∫∫ ∫∫

V V V i i i

2 2 2

1 R R U R ϕ ϕ ω χ

⊥

= ∇ ×∇ + ∇ + ∇ V

{ } ( )

{ }

( ) ( ) ( )

1 1 L p p μ γ μ = ∇× ∇× × × ⎡ ⎤ ⎣ ⎦ + ∇× × ∇× × +∇ ∇ + ∇ ⎡ ⎤ ⎣ ⎦ V V B B B V B V V i i

11 12 13 21 22 23 31 32 33

δW ( ) δW ( , ) δW ( , ) δW ( , ) δW ( ) δW ( , ) δW ( , ) δW ( , ) δW ( , ) , ,

i i i i i i i i i

U U U ν ν ν ν ν ω ω χ χ χ ω ν ν ν ν + + + + + +

these “energy terms” add to mass matrix to make a fully stable implicit system. same functions!

26

1 1/2 13 13 13 23 23 23 31 32 33 12 12 12 21 22 21 11 11 11 31 32 33 22 31 3 2 2 33 21 2 n n n v v v v v v v v v v v v v v v v v v v v v v v v v v v

S D R S S D D R R S D R S D S U D U f R S S S D D D R R R F R χ χ ψ ω ω

+ +

⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = + ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ i i i

The sparse matrix equation to be solved for the velocity variables take the form:

{ } { }

1/2 2 2 1 2 2

1 ( ) ( ) ( )

n n n

t L t L t p ρ θ δ ρ θ δ δ μ

+ +

⎧ ⎫ − = − + −∇ + ∇× × ⎨ ⎬ ⎩ ⎭ V V B B

• Corresponds to projections of the operator equation derived on earlier vg:
• Also contains 2 non-trivial sub-systems (reduced MHD) that

conserve appropriate “energy” and are numerically stable

[ ] [ ] [ ]

11 1 11 1 1 v v v n n n

S U D U R ψ

+

⎡ ⎤ ⎡ ⎤ ⎡ ⎤ = + ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ i i i

etc.

( )

2 11 11 11

et , ( ) ( , c ) .

v i i

S D U t W U

ν

ρ ν θδ δ ν = = − Sv matrix is self-adjoint!

27

First 3 δWij terms

[ ] [ ]

( )

[ ] [ ]

( )

[ ]

( )

[ ] ( )

* 4 4 4 4 2 2 2 2 11 * 2

ˆ ˆ ˆ ˆ , , 1 1 2 ˆ ˆ ˆ ( , ) , , , , , , , , , ,

i i i i iZ i i i

F F F F F F U v U v v U U v R W U F U v v U R R R R U R R R R ψ ψ ψ ψ ψ ψ δ ν ψ ν = + ⎡ ⎤ ⎡ ⎤ + − Δ − ⎣ ⎦ ⎢ ⎥ ⎣ ′ ′ ⎛ ⎞ ⎛ ⎞ ′ ′ ′ − − Δ − ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ ⎦

[ ] [ ]

( )

[ ] ( )

( )

2 4 12 * 2 2

ˆ ˆ ˆ , , , , 2 , ) , ( ˆi

i Z i i i

F v v v R R R F R W ω ω ψ ψ ψ ψ ω δ ν ψ ω ψ ν ω ′ ′ ′ ′ − Δ + − = +

( ) [ ]

( )

( ) [ ] ( )

( )

( ) [ ]

( )

* 2 * 2 2 4 2 2 3 2 2 4 1

1 ˆ ˆ ˆ ˆ , , , , , 1 1 2 ˆ ˆ ˆ , , , , , ( , ) ,

i i i i i i Z i i

F F F F v v v F v R R R F F v v R W R R R R R χ ψ ψ ψ χ ψ ν ψ χ δ ν χ χ ψ ψ χ χ

⊥

− + ′ ⎛ ⎞ ′′ ′ ′ ⎡ ⎤ − + − Δ − ⎜ ⎟ ⎣ ⎦ ⎝ ⎠ Δ + ∇ ∇ ⎡ ⎣ = ⎤ ⎦ i [ ] [ ] ( ) ( )

1 , ,

Z R R Z R R Z Z

a b a b a b a b R a b a b a b a b ϕ ≡ ∇ ×∇ ∇ = − ≡ ∇ ∇ = + i i

f f ϕ ′ ≡ ∂ ∂

• present in
• 3D

2D

• nly

note: at most second order derivatives on each scalar compatible with C1 elements

28

Magnetic Field

2

ˆ ln R f F Rz ϕ ψ ϕ = ∇ ×∇ + ∇ − A

2

1 (gauge condition) R

⊥

∇ = A i

*

f F f F ψ ϕ ϕ ψ ϕ ϕ

⊥

= ∇× ′ = ∇ ×∇ −∇ + ∇ ′ = ∇ ×∇ −∇ + ∇ B A

2 2 2

* F F R f F F R f F f

⊥

≡ + ∇ ∇ ′′ ≡ + ∇ = + i

f f ϕ ′ ≡ ∂ ∂

* * 2

1 F R ϕ ψ ψ ϕ

⊥

≡ ∇× = ∇×∇× ′ = ∇ ×∇ + ∇ − Δ ∇ J B A

2 scalar variables and a gauge condition

29

Magnetic Field Advance Equations

2 2 2

ˆ ln R f F Rz f F F F R f ϕ ψ ϕ ψ ϕ ϕ

⊥ ⊥

= ∇ ×∇ + ∇ − ′ = ∇ ×∇ − ∇ + ∇ = + ∇ A B

( )

[ ]

1 t η ∂ = ∇× × − + ⋅⋅⋅ ∂ B V B J

2 2

(1) (1)

i i

d R d R ν ϕ ν ϕ

⊥ ⊥

∇ ∇ × → ∇ ×∇

∫∫ ∫∫

i i

2 2

(1) (1)

i i

d R d R ν ϕ ν ϕ ∇ → ∇

∫∫ ∫∫

i i

2 2

(1) (1)

i i

d R d R ν ν

⊥

− ∇ → ∇

∫∫ ∫∫

i i

Not needed!

Energy Conservation:

( ) ( ) ( ) ( ) ( ) [ ]

1 1 ( ) 1 1 1 p t t p p p p t p p γ γ γ γ γ ∇× × ⎡ ⎤ ⎣ ⎦ − ∇× ∂ ⎡ ⎤ = ∇× × − ∇ = − ⎢ ⎥ ∂ ⎣ ⎦ ∂ ⎡ ⎤ = ∇× × = + ∇ × × ⎡ ⎤ ⎣ ⎦ ⎢ ⎥ ∂ ⎣ ⎦ ∂ = − ∇ × ⎡ ⎤ ⎣ + ∇ = − ∇ + − ∂ − ∇ ⎦ − ∇ V V B B B B V B V V B B V V V B V V B V B B i i i i i i i i i i

For energy conservation, the like colored terms must cancel exactly. Since this only requires that the projections we take of the momentum equation are equivalent to the dot product with the velocity, we will have energy conservation for each of the 3 velocity fields:

2 2 2 2 2 2

1 U U R U R R R R R ϕ ϕ ω ϕ ϕ ω ϕ χ

⊥

= ∇ ×∇ = ∇ ×∇ + ∇ = ∇ ×∇ + ∇ + ∇ V V V

Reduced MHD Full MHD

31

2-Variable 3D Toroidal subset of full equations….or, (1,1) component ( )

[ ]

* * * * 2 2 2 4 2 2 2 4 * * * 2 2 4 2 2

1 1 1 1 , , , 1 1 1 1 , ( ) F F U U U U U R R R R R R R R F U U R R R R R ψ ψ ψ ψ μ μ ψ ψ η ψ η ψ

∗ ∗ ⊥ ⊥

⎡ ⎤ ⎡ ⎤ Δ Δ ⎛ ⎞ ⎡ ⎤ ′ ′ ′′ Δ = − + + Δ + + Δ Δ + ∇ ∇ ⎜ ⎟ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎝ ⎠ ⎣ ⎦ ⎣ ⎦ ⎡ ⎤ ′ ′′ Δ = − Δ + ∇ ∇ + ∇ ∇ Δ + ∇ ⎢ ⎥ ⎣ ⎦

• i
• i

i

2

ˆ ln F RZ R U ψ ϕ ϕ = ∇ − = ∇ ×∇ A V

( ) ( )

[ ] [ ] ( )

* * * 2 2 * 2 2 2 4 2 4

1 1 , ( ) , , , ,

D i i i i i i i

F U U U t W U U U R R R R R R ψ ν θδ δ ν ν ν ψ ν ψ μ ν ν ⎡ ⎤ Δ Δ Δ ′ ′′ − − = − − + Δ − ∇ ∇ ⎢ ⎥ ⎣ ⎦

• i

( ) [ ] [ ]

( )

[ ]

( )

[ ]

( )

( ) [ ] [ ]

2 4 4 2 2 * * 6 2 4

1 , , , , , , , , , , , ,

D

F F U V U V V U R R R W U V RdRdZ F F U V V U V U R R R ψ ψ ψ ψ δ ψ ψ ψ ⎧ ⎫ ′ ′ + − ⎪ ⎪ ⎪ ⎪ = ⎨ ⎬ Δ ⎪ ⎪ ′′ ′ ⎡ ⎤ − − − Δ ⎣ ⎦ ⎪ ⎪ ⎩ ⎭

∫∫

These reduced equations have an energy conservation theorem, and the ideal terms have an associated energy principle and variational form Weak form for the velocity equation:

32

4-Variable 3D Toroidal subset of full equations….or, 2x2 submatrix

2 2 2

ˆ ln R f F Rz R U R ϕ ψ ϕ ϕ ω ϕ = ∇ ×∇ + ∇ − = ∇ ×∇ + ∇ A V

( )

( ) ( ) ( ) ( )

2 2 2 2 2 2 2 2 2 2 2 2 4

, , , , , , U U R U U U R R f f R z R F F F F R f U R U R R R z R R ψ ω ω ω ψ ψ ψ μ ψ ψ μ

∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗

⎡ ⎤ ⎡ ⎤ Δ ∂ Δ ′ ′ ′ ′ Δ = − − − Δ − + − Δ Δ − Δ ⎢ ⎥ ⎢ ⎥ ∂ ⎣ ⎦ ⎣ ⎦ ⎡ ⎤ ∂ ⎛ ⎞ ⎡ ⎤ ⎛ ⎞ ⎡ ⎤ ′ ′ ′′ ′′ + Δ + + + + Δ Δ + ∇ ∇ ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ ⎢ ⎥ ⎢ ⎥ ∂ ⎝ ⎠ ⎣ ⎦ ⎝ ⎠ ⎣ ⎦ ⎢ ⎥ ⎣ ⎦

• i

( ) [ ] (

)

2 2 2 * * 2 C 2

1 , , , , , 2 R R U R F f f F p R R ω ω ωω ψ ψ ψ ψ μ ω μ ω

∗

′ ′ ′ ′ ′ ′ ′′ ⎡ ⎤ ⎡ ⎤ = − − − + + − − + Δ + ⎣ ⎦ ⎣ ⎦

• [

] ( )

( )

( ) [ ]

2 * * * 2 2 2 2 2 * * * * 2 2 2 2

1 , , , 1 1 1 1 1 , , , F U R f F U U f F R R F F R R R F f f F f R R R R ω ψ ω ϕ η ψ η ϕ ψ η ψ ψ ψ ψ ψ ψ ψ ψ ψ ψ ϕ

⊥ ⊥ ⊥ ⊥ ⊥ ⊥ ⊥ ⊥ ⊥ ⊥ ⊥

′ ⎡ ⎤ ′ ′ ⎡ ⎤ ∇ − ∇ − ∇ ×∇ + ∇ + ∇ ×∇ ′ ⎡ ⎤ − − + Δ + ⎢ ⎥ ⎣ ⎦ ⎢ ⎥ ⎢ ⎥ ∇ ∇ = ∇ ∇ + ∇ ⎢ ⎥ ⎡ ⎤ ⎢ ⎥ ′ ′ ′ ′ + + + ′ ′ − Δ ∇ − ∇ + ∇ − Δ ∇ ×∇ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦

• i

i i

* 2 2 2 2 2 2 * * * 2 2 4 2

1 1 1 F U f F R R R F R f R F F F f R R R R η ϕ ω ϕ ψ ω η ϕ ψ ψ ϕ ψ ϕ ψ ψ

⊥ ⊥ ⊥ ⊥ ⊥ ⊥ ⊥ ⊥ ⊥ ⊥

⎡ ⎤ ′ ′ ∇ ×∇ − ∇ ×∇ − ∇ + ∇ + ∇ ×∇ ⎢ ⎥ = ∇ = ∇ ⎢ ⎥ ⎢ ⎥ ′ ′ − Δ ∇ ×∇ − ∇ ×∇ + ∇ − Δ ∇ ⎢ ⎥ ⎣ ⎦

• i

These equations also have an energy theorem and associated energy principle.

33

3D C1 elements by combining Q18 triangles in (R,Z) Hermite Cubic representation in the toroidal angle φ

( ) ( ) ( )

2 2 1 2

( ) | | 1 2 | | 1 ; ( ) | | 1 x x x x x x Φ = − + Φ = −

]

18 1 2 , 1 , 2 1 1 2 , 1 1 , 1 2

( , , ) ( , ) ( / ) ( / ) ( / 1) ( / 1)

j j k j k j j k j k

U R Z R Z U h U h U h U h ϕ ν ϕ ϕ ϕ ϕ

= + +

⎡ = Φ + Φ ⎣ + Φ − + Φ −

∑

Each toroidal plane has two Hermite cubic functions associated with it Solution for each scalar function is represented in each triangular wedge as the product of Q18 and Hermite functions All DOF are still located at nodes: => very efficient representation

34

3D Nonlinear Solver Strategy

• In 2D, solve efficiently with direct solver up to (200)2 nodes
• In 3D, leads to block triangular structure

1 1 1 1 1 1 1 1 1

. . . . . .

j j j j j j N N

x y x y x y x y x y

− − + +

⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦

j j j N N N

B C A A B C C A B i i i i i i i i i i i i i i i i Block Jacobi preconditioner corresponds to multiplying each row by PETSc has the capability of doing this using SuperLU_Dist

• 1

j

B

, ,

j j j

A B C

are 2D sparse matrices at plane j

N=1 Resistive Internal Kink mode in CMOD with S=107

Perturbed Current with Mesh Close-up

S=108 S=107 S=106 S=105

• M3D-C1, is now being used for linear physics

studies in NSTX, CMOD and ITER

• high order C1 finite elements, adaptive mesh,

and fully implicit time advance allow high resolution studies of localized modes

• Now being used to study tearing (and double

tearing) modes at realistic S values, including pressure (Glasser) stabilization

(Top) Equilibrium current density with adaptive mesh superimposed. (Left) perturbed current density for (1,1) tearing mode at different S. Rightmost figure corresponds to NSTX parameters

36

High-S tearing mode studies

q0=0.85 q0=0.60 η=10-5 η=10-7 η=10-5 η=10-7

• Poloidal Vorticity Δ*U
• Toroidal Current Density Δ*ψ
• CMOD:

β0=.006

n=1 Double Tearing Mode in NSTX….S = 108

q-profile Toroidal current Vorticity Normal displacement

39

40

Linear ELM1s: Code Verification

We have performed detailed benchmarking for ELM unstable equilibrium in the ideal limit between M3D-C1 and GATO and ELITE up to n=40

• required discontinuous η and ρ profiles

with jump of 108

Ferraro

ideal limit: plasma resistivity 0. vacuum region resistivity ∞, vacuum region density 0

1Edge Localized Modes

Studies have been extended to:

• diverted equilibrium (JT-60)
• finite resistivity in the plasma and SOL
• realistic density profiles

Close-up showing M3D-C1 triangular adaptive mesh

Ferraro 41

Linear ELMs: Code Verification-2

Comparison of eigenfunctions of normal plasma displacement for “ideal limit” and more realistic Spitzer resistivity with SOL with M3D-C1

Ideal MHD limit Sptizer resistivity with SOL

42

Linear ELMs: Code Verification-3

Plasma burst outboard, midplane n reduced Density to outboard divertor Inboard edge instability Density to inboard divertor γ

Outer div Inner div time (τA)

Multi-stage ELM – DIII-D 119690

Poloidal rotation (?)

Sugiyama

Linear mode growth and mode consolidation

43

Non-linear ELM simulation with M3D

Full simulation of nonlinear ELM event shows complex structure with secondary instabilities

Initially many unstable linear modes. These rapidly consolidate into lower-n field-aligned mode ``filaments'' (n=6-10 at t=43) Similar to what is seen experimentally. n pert n ψ-pert u RJφ

First 50 τA: linear mode growth Nonlinear harmonic consolidation

Sugiyama 44

Early time: T and n ballooning in rapid burst

T n t=21.5 τA 42.8 62.3 83.4 104.6

Sugiyama 45

Longer time: T t=43 126 227 461 529

Sugiyama 46

Longer time: n t=43 126 227 461 529

Sugiyama 47

ECCD Stabilization of NTM

NIMROD code calculates the MHD growth of NTM c GENRAY code computes wave induced ECCD current drive term

Code coupling provided by SWIM framework

Jenkins 48

Results to date are for an equilibrium that is tearing unstable and using a model toroidally localized CD term

Close up Model current drive source applied to original O-point in 1 toroidal location. (2,1) island shrinks, becomes (4,2) (4,2) island shrinks, (2,1) grows New (2,1) 900 our of phase with old

Jenkins 49

(2,1) islands grow up again after (4,2) island has been suppressed; new islands are 90° out

• f phase (poloidally)

from the old ones.

As RF suppresses the original islands, new islands arise

Jenkins 50

Study of saturated mode in NSTX-Motivation

NSTX shot 124379 has a steadily growing 2,1 mode with no apparent trigger seen by the USXR, Dα, or neutron diagnostics.

Gerhardt

51

Eigenfunction Analysis of Multichord Data Suggests Coupling to 1,1 Ideal Kink

2,1 only 2,1 + 1,1 pert

Gerhardt

52

M3D simulation of saturated mode in NSTX when q0 > 1

Saturated n=1 mode can set develop when q0 slightly > 1, as seen in Poincare plot on left. Can flatten temperature (right) and also drive m=2 islands. Breslau, et al. IAEA 2010

54

VDE1 and Plasma Disruption simulations in ITER

1Vertical Displacement Event

(a) Poloidal flux, (b) toroidal current, and (c) temperature during a vertical displacement event. A VDE brings the plasma to the upper wall where a (m,n) = (1,1) kink mode grows. Forces on the vacuum vessel are calculated.

55

Runaway electron evolution in disrupting plasma is computed.

Simulation of DIII-D Ar pellet experiments. Runaway electrons of different energy shown. Synchrotron emission on right.

Izzo

Thin Resistive Wall Boundary Conditions

R Z

ϕ

ˆ n ˆ t

P V Thin resistive wall Tangential integration path R

ϕ P V V (a) (b)

New instabilities can be present if the plasma is surrounded by a thin resistive wall. This can be modeled by modifying the boundary conditions. All the boundary conditions follow from imposing that the normal component of B is continuous across the wall, and the tangential components can have a jump. Follows from: ∇ = B i

Thin Resistive Wall Boundary Conditions-2

ˆ ˆ

V P

n n ⋅ = ⋅ B B

2 V V V

φ φ = ∇ ∇ = B

[ ]

1 ˆ ˆ ˆ ˆ ˆ

W W P V P V

n l l t l R η η ϕ ϕ δ ϕ δ ∂ ∂ ∂ ⎡ ⎤ ⋅ = ⋅ − ⋅ + ⋅ − ⋅ ⎣ ⎦ ∂ ∂ ∂ B B B B B magnetic field on vacuum side of wall magnetic field on plasma side of wall magnetic scalar potential in wall resistivity of wall thickness of wall

V P V W

φ η δ = = = = = B B

Without the presence of flow in the plasma, a wall stabilized plasma will become unstable as wall resistance is increased

With sufficient plasma flow, the unstable mode will stabilize due to doppler shifted resonance with sound continua

• Implicit Methods

– Use differential approximation to reduce matrix size and improve condition number

• Highly magnetized plasma

– Stream-function/potential representation of velocity and magnetic fields

• Momentum equation projections

– Gives energy conservation for full and reduced equation sets

• Finite Elements

– High-order C1 elements needed for high-order equations.

• Solver Strategy

– Block Jacobi preconditioner

• Recent Results

– Demanding ELM ideal MHD benchmarking studies give excellent results – Lundquist numbers up to 108 or higher are possible

60

Summary

61