The Galerkin Finite Element Method Stephen C. Jardin Princeton - - PowerPoint PPT Presentation

MHD Simulations for Fusion Applications

Lecture 3

The Galerkin Finite Element Method

Stephen C. Jardin Princeton Plasma Physics Laboratory

CEMRACS ‘10 Marseille, France July 21, 2010

1

Outline of Remainder of Lectures

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

differencing

Today Tomorrow

• 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

Weak Form of a Differential Equation

Consider the ordinary differential equation on the domain [a,b]: ( ) ( ) ( ) ( ) (1) d d p x q x U x f x dx dx ⎡ ⎤ ⎛ ⎞ − + = ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎣ ⎦ We choose a function space H1, and multiply (1) by an arbitrary member of that function space, and integrate over the domain: ( ) ( ) ( ) ( ) ( ) ( ) (2 )

b b a a

d d V x p x q x U x dx V x f x a dx dx ⎡ ⎤ ⎛ ⎞ − + = ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎣ ⎦

∫ ∫

( ) ( ) ( ) ( ) (2 )

b b a a

dV dU p x q x VU dx V x f x b dx dx ⎡ ⎤ ⎛ ⎞ + = ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎣ ⎦

∫ ∫

Or, after integrating by parts (assuming the boundary term vanishes): If we require that Eq. (2) hold for every function V(x) in the function space H1, it is called the weak form of Eq. (1). Note that only first derivatives occur in (2b), whereas 2nd derivatives occur in (1)

Galerkin Finite Element Method

The Galerkin finite element method is a discretization of the weak form

• f the equation. To apply, we chose a finite dimensional subspace S of

the infinite dimensional Hilbert space H1, and require that Eq. (2) hold for every function in S.

1

( ) ( ) ( ) ( ) 1,

b b N j i j i j i j a a

d d a p x q x dx x f x i N dx dx ϕ ϕ ϕ ϕ ϕ

=

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

∑ ∫ ∫

1

( ) ( )

N j j j

U x a x ϕ

=

=∑

( ) ( ) ( ) ( )

b j i ij i j a b i i a ij j i

M a d d M p x q x dx dx dx b x f x b dx ϕ ϕ ϕ ϕ ϕ ⎡ ⎤ ⎛ ⎞ = + ⎢ ⎥ ⎜ ⎟ ⎝ ⎠ ⎣ ⎦ = =

∫ ∫

known basis functions that span S aj are amplitudes to be solved for Letting Equation (2b) becomes Or, In the finite element method, the basis functions are low order polynomials ( ) ( )

i

V x x ϕ =

Simple Example: Linear Elements in 1D

1

( ) ( )

N j j j

U x a x ϕ

=

= ∑

( )

i x

ϕ

1( ) i

x ϕ +

1( ) i

x ϕ −

xi xi-1 xi+1 L 1

, (constants) ( ) p p q q f f x = = =

, at i i j j

x x ϕ δ = =

linear: xi = ih h=L/N

,

2 1 1 2 1 1 1 2 1 1 2

i j i j

p dx K p h ϕ ϕ − ⎡ ⎤ ⎢ ⎥ − − ⎢ ⎥ ∇ ∇ = = ⎢ ⎥ − − ⎢ ⎥ − ⎣ ⎦

∫∫

i

,

4 1 1 4 1 1 4 1 6 1 4

i j i j

h q v v dx A q ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ = = ⎢ ⎥ ⎢ ⎥ ⎣ ⎦

∫∫

“mass matrix” “stiffness matrix”

Leads to sparse, diagonally dominant matrices: Mi,j = Ki,j + Ai,j

i i

b fdx ϕ = ∫∫

ij j i

M a b =

1

( ) ( ) ( ) ( ) 1,

b b N j i j i j i j a a

d d a p x q x dx x f x i N dx dx ϕ ϕ ϕ ϕ ϕ

=

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

∑ ∫ ∫

Linear elements only

• verlap with

their nearest neighbors

Note: We were able to solve a 2nd order ODE with linear elements!

1

( ) ( )

N j j j

U x a x ϕ

=

= ∑

( )

i x

ϕ

1( ) i

x ϕ +

1( ) i

x ϕ −

xi xi-1 xi+1 L 1

1

( ) ( ) ( ) ( ) 1,

b b N j i j i j i j a a

d d a p x q x dx x f x i N dx dx ϕ ϕ ϕ ϕ ϕ

=

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

∑ ∫ ∫

( ) ( ) ( ) ( ) d d p x q x U x f x dx dx ⎡ ⎤ ⎛ ⎞ − + = ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎣ ⎦ ( )

b i a

dx x ϕ

∫

ij j i

M a b =

Extension to Higher Order Equations

Suppose we have a 4th order equation:

4 4

( ) ( ) d U x f x dx = Form the weak form:

4 4 3 3 2 2 2 2

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

b b i i a a b b i i a a b b i i a a

d x U x dx x f x dx dx d d U dx x f x dx dx dx d d U dx x f x dx dx dx ϕ ϕ ϕ ϕ ϕ ϕ = − = =

∫ ∫ ∫ ∫ ∫ ∫

(assumes the boundary terms vanish) Application of the Galerkin method to a 4th order equation requires that the second derivative of the basis functions be defined everywhere. For a 2nd order equation, only the first derivative need be defined.

1

( ) ( )

N j j j

U x a x ϕ

=

= ∑

Integrate by parts twice:

Elements with higher continuity

Linear elements are continuous, but only the function and first derivative are well defined. To represent a second derivative, both the function and first derivative must be continuous. They must be a subset of H2 j-1 j j+1 j-1 j j+1

1 ( ) j x

ϕ

2 ( ) j x

ϕ

The Hermite Cubic elements associate two basis functions with each node j, defined on the interval xj-1 < x < xj+1 as:

( ) ( ) ( )

2 1 2 2 1, 1 2, 2

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

j j j j

y y y y y y x x x h x x x h h ϕ ϕ ϕ ϕ ϕ ϕ = − + = − − ⎛ ⎞ = ⎜ ⎟ ⎝ ⎠ − ⎛ ⎞ = ⎜ ⎟ ⎝ ⎠

These have the property that:

• 1. In the interval xj < x < xj+1 , a complete

cubic can be represented by:

• 2. The first derivatives vanish at the

element boundaries 3. The second derivatives are defined everywhere

1, 2, 1, 1 2, 1

( ), ( ), ( ), ( ),

j j j j

x x x x ϕ ϕ ϕ ϕ

+ +

j-1 j j+1 j-1 j j+1

1 ( ) j x

ϕ

2 ( ) j x

ϕ

Hermite Cubic Finite Elements

1, 2, 1

( ) ( ) ( )

N j j j j j

U x a x b x ϕ ϕ

=

⎡ ⎤ = + ⎣ ⎦

∑

In the interval [ j , j+1 ] , the function is written as:

( )(

)

( )(

)

1, 2, 1 1, 1 1 2, 1 2 1 1 3 1 1 2 3 1 2 3

( ) ( ) ( ) ( ) ( ) 3 2 3 2 2

j j j j j j j j j j j j j j j j j j

U x a x b x a x b x a b x a hb a hb x h a hb a hb x h c c x c x c x ϕ ϕ ϕ ϕ

+ + + + + + + +

= + + + = + + − − + − + + − + = + + +

Simple physical interpretation: aj is value of function at node j bj is value of derivative at node j

1 2 2 1 2 3 2 3 2 1 3

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

j j j j

a c b c a c h h h h b c h h h h

+ +

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

Complete cubic polynomial determined by values of function and derivative at endpoints

j j+1

1 ( ) j x

ϕ

2 ( ) j x

ϕ

Hermite Cubic Finite Elements-2

1, 2, 1

( ) ( ) ( )

N j j j j j

U x a x b x ϕ ϕ

=

⎡ ⎤ = + ⎣ ⎦

∑

In the interval [ j , j+1 ] , the function is written as:

( )(

)

( )(

)

1, 2, 1 1, 1 1 2, 1 2 1 1 3 1 1 2 3 1 2 3

( ) ( ) ( ) ( ) ( ) 3 2 3 2 2

j j j j j j j j j j j j j j j j j j

U x a x b x a x b x a b x a hb a hb x h a hb a hb x h c c x c x c x ϕ ϕ ϕ ϕ

+ + + + + + + +

= + + + = + + − − + − + + − + = + + +

1 1( ) j

x ϕ

+ 2 1( ) j

x ϕ

+

Complete cubic in each interval!

Linear elements 1st Hermite cubic 2nd Hermite cubic

( ) x ϕ ( ) x ϕ′ ( ) x ϕ′′

j-1 j j+1 j-1 j j+1 j-1 j j+1 not defined

Comparison of elements and derivatives

Types of Finite Elements in 2D

Shape Order Continuity

,

( , )

i k ik i k M i k

x y C x y φ

+ < ≥

= ∑

rectangle triangle Order M if complete polynomial to that

• rder

C0: continuous

across element boundaries

C1: continuous

first derivatives across element boundaries

DG: discontinuous

Galerkin

Element Shape

Rectangle:

• natural for structured mesh
• allows tensor product forms

bi-linear: t=1, bi-quadratic: t=2, etc

• cannot refine locally without

hanging nodes Triangle:

• natural for unstructured mesh
• easier to fit complex shapes
• can easily refine locally

,

( , )

i k ik i k t

x y a x y φ

≤ ≤

= ∑

Element Order

If an element with typical size h contains a complete polynomial of order M, then the error will be at most of order hM+1

This follows directly from a local Taylor series expansion: Thus, linear elements will be O(h2) quadratic elements will be O(h3) cubic elements will be O(h4) quartic elements will be O(h5) complete quintic elements will be O(h6)

1 ,

1 ( , ) ( ) ( ) ( ) !( )!

k M k l k l M l k l k l x y

x y x x y y O h l k l x y φ φ

− + − = =

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

∑∑

Element Continuity

Theorem: A finite element with continuity Ck-1 belongs to Hilbert space Hk, and hence can be used for differential operators with order up to 2k Continuity Hilbert Space Applicability C0 H1 second order equations C1 H2 fourth order equations

This applicability is made possible by performing integration by parts in the Galerkin method, shifting derivatives from the unknown to the trial function

[ ]

( )

2 2 2 2

recall: ( , ) ( , ) ( , ) ( , )

i i domain domain i i domain domain

v f x y dxdy f x y v dxdy v f x y dxdy f x y v dxdy φ φ φ φ ∇ ∇ = − ∇ ∇ ⎡ ⎤ ∇ ∇ = ∇ ∇ ⎣ ⎦

∫∫ ∫∫ ∫∫ ∫∫

i i

Hk means that derivatives exist up to order k

The vast majority of the literature concerns C0 elements,. Here we concentrate on C1 elements because we are interested in higher order equations

NOTE: requires the trial function have appropriate boundary conditions

Reduced Quintic 2D Triangular Finite Element: Q18

θ P1(x1,y1) P3(x3,y3) P2(x2,y2)

• rigin

ξ η

20 1

( , )

k k

m n k k

a φ ξ η ξ η

=

= ∑

k mk nk 1 2 1 3 1 4 2 5 1 1 6 2 7 3 8 2 1 9 1 2 10 3 11 4 12 3 1 13 2 2 14 1 3 15 4 16 5 17 3 2 18 2 3 19 1 4 20 5 For C1, require that the normal slope along the edges φn have only cubic variation: 5b4ca16 + (3b2c3 – 2b4c)a17 + (2bc4 – 3b3c2) a18 + (c5 – 4b2c3)a19 – 5bc4a20 = 0 5a4ca16 + (3a2c3 – 2a4c)a17 + (-2ac4 – 3a3c2) a18 + (c5 – 4a2c3)a19 – 5ac4a20 = 0 20 – 2 = 18 unknowns:

These are determined in terms of [ φ, φx, φy , φxx, φxy, φyy ] at P1,P2,P3

Implies C1 continuity at edges and C2 at nodes ! x y

Specify the function and it’s first and second derivatives at the 3 nodes (φ, φx, φy, φxx, φxy, φyy) Guarantees C0 continuity Note: general quintic has 21 terms. 1 has been set to zero, and 2 others will be determined by constraints.

ξ,η are local

• rthogonal

coordinates

No mk=4, nk=1 term

6 Unknowns at each node!

2 3 4 5 1 2 3 4 1 2 3 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3

1 1 2 3 4 5 1 b b b b b b b b b b b b

ξ η ξξ ξη ηη ξ η ξξ ξη ηη ξ η ξξ ξη ηη

φ φ φ φ φ φ φ φ φ φ φ φ φ φ φ φ φ φ − − − ⎡ ⎤ ⎢ ⎥ − − ⎢ ⎥ − − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦

2 3 2 2 3 2 3 4 5 2 3 4 2 3 2 3 2

2 6 12 20 1 2 3 2 2 2 2 1 1 2 3 4 5 1 2 6 12 20 1 2 3 b b b b b b b b a a a a a a a a a a a a a a a a a − − − − −

2 3 2 3 4 5 2 3 4 2 3 4 2 3 2 3 2 3 4

2 2 2 2 1 1 1 2 3 4 5 2 2 2 2 1 2 3 4 2 6 12 20 5 a a a c c c c c c c c c c c c c c c c c c c c c c a

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 2 3 4 5 4 18 4 3 2 2 3 19 2 3 4 5 4 4 20 4 3 2 2 3

3 2 5 2 3 4 3 2 5 5 2 3 4 a a a a a a a a a a a a a a a a a a c ac c a c ac a c a c a c a b c bc c a b c bc b c b c b c ⎡ ⎤ ⎡ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − − + − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − − − − ⎣ ⎢ ⎥ ⎣ ⎦ ⎤ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎦

[ ] [ ][ ]

T A ′ Φ =

• r,

20 1

( , )

k k

m n k k

a φ ξ η ξ η

=

=∑

using the expansion Calculate the value of φ and it’s derivatives wrt local coords. ξ and η at the 3 vertex points Combine with the 2 constraint equations:

[ ] [ ]

1

A T − ′ ⎡ ⎤ = Φ ⎣ ⎦ How to calculate the coefficients in terms of the unknowns:

How to calculate the coefficients in terms of the unknowns-2:

The previous viewgraph calculated the coefficient matrix in terms of derivatives in the rotated coordinate system, which is different for each

• element. To get the coefficient matrix in terms of the vector containing the

actual derivatives with respect the Cartesian coordinates (x,y), we have to apply a rotation matrix that depends on the triangle’s orientation.

⎡ ⎤ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎣ ⎦

1 1 1

R R R R

2 2 2 2 2 2

1 cos sin sin cos cos 2sin cos sin sin cos cos sin sin cos sin 2sin cos cos θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ θ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − − ⎢ ⎥ − ⎣ ⎦

1

R 2

G = T R i

Relates the coefficient matrix directly to the function and derivatives wrt (x,y) at the nodes 20x18 matrix consisting of the first 18 columns of T-1 Each triangle has it's own 18 18 gij × G

20 1

i i

m n j ij i

v g ξ η

=

=∑

The Trial Functions:

20 20 18 18 1 1 1 1

i i i i

m n m n i ij j j j i i j j

a g v φ ξ η ξ η

= = = =

= = Φ = Φ

∑ ∑∑ ∑

These are the trial

• functions. There are

18 for each triangle. The 6 shown here correspond to one node, and vanish at the other nodes, along with their derivatives Each of the six has value 1 for the function or one of it’s derivatives at the node, zero for the

• thers.

Note that the function and it’s derivatives (through 2nd) play the role of the amplitudes: 6 at each node!

ai = gij Φj Function and derivatives wrt (x,y)

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

( , )

k k k k k l k l k l k l

m n m n k k k k m m n n m m n n k l kl kl k l k l k l

a m n a a K K m m n n φ ξ η ξ η ξ ξ η η φ ξ η ξ η

− − = + − + + + − = =

⎡ ⎤ ∇ = ∇ + ∇ ⎣ ⎦ ∇ = = +

∑ ∑∑

1 1 1

( ) ! ! ( , ) ( 2)! ( 2, ) ( , 2)

m m m n n element kl kl k l k l k l k l k l k l element

a b m n d d c F m n m n K K d d m m F m m n n n n F m m n n ξ η ξ η ξ η

+ + + ⎡

⎤ − − ⎣ ⎦ = ≡ + + ≡ = + − + + + + −

∫∫ ∫∫

Integrals involving powers of the local coordinates can be done analytically: Thus, for example: ai = gijΦj

All matrix elements easily calculated in terms of the ai (or gij)

20 1

i i

m n j ij i

v g ξ η

=

=∑

20 20 18 18 1 1 1 1

i i i i

m n m n i ij j j j i i j j

a g v φ ξ η ξ η

= = = =

= = Φ = Φ

∑ ∑∑ ∑

Very compact representation compared to a popular C0 Element because all unknowns on nodes

6 6 6 Lagrange Cubic: C0, h4 Reduced Quintic: C1, h5 6 6 6 6 6 new unknowns: 2 new triangles 6/2 = 3 unknowns/ triangle 9 new unknowns: 2 new triangles 9/2 = 41/2 unknowns/ triangle split split

Comparison of reduced quintic to other popular triangular elements

Vertex nodes Line nodes Interior nodes accuracy

• rder hp

UK/T continuity linear element 3 2 ½ C0 Lagrange quadratic 3 3 3 2 C0 Lagrange cubic 3 6 1 4 4½ C0 Lagrange quartic 3 9 3 5 8 C0 reduced quintic 18 5 3 C1 *

N M Linear Elements Reduced Quintic Elements

2

1

L

E a N =

5

1

Q

E b M =

15 2 / 5 2 / 5

same ~ error

L Q

b E E M N N a ⎛ ⎞ => = => = ⎜ ⎟ ⎝ ⎠

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . x x x x x x x x x x x x x x x x x x x x x ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦

. . . . . . . . x x x x x x x x x x x x x ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦

N2 36N4/5 2N 12N2/5 Win for N > 20 !

2

S ∇ Φ =

Second order equation:

N M Linear Elements Reduced Quintic Elements

2

1

L

E a N =

5

1

Q

E b M =

15 2 / 5 2 / 5

same ~ error

L Q

b E E M N N a ⎛ ⎞ => = => = ⎜ ⎟ ⎝ ⎠

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . x x x x x x x x x x x x x x x x x x x x x ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦

. . . . . . . . x x x x x x x x x x x x x ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦

4N2 36N4/5 4N 12N2/5 Win for N > 6 !

4

S ∇ Φ =

2

∇ Φ = Ψ

2

S ∇ Ψ =

Fourth order equation:

Anisotropic Diffusion

2

BB S t B φ κ φ κ φ ∂ ⎡ ⎤ = ∇ ∇ + ∇ ∇ + ⎢ ⎥ ∂ ⎣ ⎦

• i

i i

18 1 j j j

v φ

=

= Φ

∑

20 ( ) ( ) 1

i i

m n l l j ij i

v g ξ η

=

=∑

Apply Galerkin method: Typically κ|| >> κ⊥

( ) ( ) ( ) ( ) 2 ( ) ( ) ( ) 2 l l l l j j j j l l l j j j

BB v dxdz v v v S d d t B v BB v v S d d B φ κ φ κ φ ξ η κ φ κ φ ξ η ⎡ ⎤ ⎡ ⎤ ∂ = ∇ ∇ + ∇ ∇ + ⎢ ⎥ ⎢ ⎥ ∂ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎧ ⎫ = − ∇ ∇ − ∇ ∇ + ⎨ ⎬ ⎩ ⎭

∫∫ ∫∫ ∫∫

• i

i i i i i

Now,

( ) ( ) ( )

, ˆ

l l l j j j

v B v z v ψ ψ ⎡ ⎤ ∇ = ∇ ×∇ = ⎣ ⎦

• i

i

[ ]

18 , 1

( , ) , ,

j j j k k k

v d d v d d R ξ η φ ξ η ψ φ ψ ξ η

=

⎡ ⎤ ∇ ∇ = ⎣ ⎦ = Φ

∫∫ ∫∫ ∑

BB i i

( )(

)

18 18 20 20 20 20 , , , , , 1 1 1 1 1 1 18 18 , , , 1 1

( 2, 2)

j k p j q i r k s l p q q p r s s r i l p q r s p q r s p q r s i l j k i l i l i l

R g g g g m n m n m n m n F m m m m n n n n P

= = = = = = = =

≡ − − × + + + − + + + − Ψ Ψ = Ψ Ψ

∑∑∑∑∑∑ ∑∑

Evaluation of the Matrix Elements

( )(

)

20 20 20 20 , , , , , , , 1 1 1 1

( 2, 2)

j k i l p j q i r k s l p q q p r s s r p q r s p q r s p q r s

P g g g g m n m n m n m n F m m m m n n n n

= = = =

≡ − − × + + + − + + + −

∑∑∑∑

[ ]

ˆ ,

x y y x

z a b a b a b a b a b

ξ η η ξ

ψ = ×∇ = − = − B

Anisotropic Diffusion

N..number of points per side

10 20 30 40 60

RMS error in steady-state

1e-6 1e-5 1e-4 1e-3 1e-2 1e-1

κ|| = 10**4 κ|| = 10**5 κ|| = 10**6 κ|| = 10**7 κ|| = 10**8

power law

N-4 N-5

2

BB S t B φ κ φ κ φ ∂ ⎡ ⎤ = ∇ ∇ + ∇ ∇ + ⎢ ⎥ ∂ ⎣ ⎦

• i

i i

Solve to steady state cos cos x y S L L π π ψ = =

Shows greater than N-5 convergence

2D Incompressible MHD

[ ]

2 2 2 4 2

, , , t t φ φ φ ψ ψ μ φ ψ ψ φ η ψ ∂ ⎡ ⎤ ⎡ ⎤ ∇ + ∇ − ∇ = ∇ ⎣ ⎦ ⎣ ⎦ ∂ ∂ + = ∇ ∂

2 2 2 2 2 4 4 2 2 1 1

, , , ,

n n n n n n n n n

t t t t t t t t t t φ φ θδ φ φ θδ φ ψ θδ ψ ψ θδ ψ μ φ θδ φ ψ ψ θδ ψ φ θδ φ η ψ θδ ψ φ φ ψ ψ φ ψ δ δ

+ +

⎡ ⎤ ⎡ ⎤ ∇ + ∇ + ∇ + − ∇ + ∇ + ⎣ ⎦ ⎣ ⎦ ⎡ ⎤ = ∇ + ∇ ⎣ ⎦ ⎡ ⎤ ⎡ ⎤ + + + = ∇ + ∇ ⎣ ⎦ ⎣ ⎦ − − = =

• 18

1 n n j j j

v φ

=

= Φ

∑

18 1 n n j j j

v ψ

=

= Ψ

∑

20 1

i i

m n j ij i

v g ξ η

=

=∑

Multiply equations by each trial function and integrate over space

“reduced MHD” φ is stream function ψ is poloidal flux

θ-centering….time centered about n+1/2 for θ=0.5

note:

Leads to the Matrix Implicit System

11 12 * * , , , , , , 21 22 * * , , , , , , * , , , , , , , , 11 12 , 21 22

[ ] [ ] [ ( (1 ) ]

j j i j i j k k i j i j k k j j i j k k i j i k j k i j n i j i j k k i j k k n i j k k i i j j j j j

S S A t G B tG S S tK M t K A A t G G t G G B D D D D θδ μ θδ θδ θδ η δ θ δ θ θ μ ⎡ ⎤ ⎡ ⎤ + Φ + − Ψ = ⎢ ⎥ ⎢ ⎥ Ψ + Φ − ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎧ ⎫ − Φ − Φ ⎪ ⎪ Ψ − ⎨ ⎬ + − ⎡ ⎤ ⎪ ⎪ ⎩ ⎭ = ⎢ ⎥ ⎢ ⎥ ⎣ ⎦

* , * 1 , , , 2 * 1 , , 2 ,

) [ ( ) ( ) (1 ) ]

j k k n i j i k j k k n i j k k k i j

M t K tK A δ θ δ θ θ η ⎡ ⎤ Ψ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎧ ⎫ − Φ − Φ ⎢ ⎥ ⎪ ⎪ − Ψ + Ψ ⎨ ⎬ ⎢ ⎥ − − ⎪ ⎪ ⎢ ⎥ ⎩ ⎭ ⎣ ⎦

11 12 1 11 12 21 22 1 21 22 n n j j j j j j n n j j j j j j

S S D D S S D D

+ +

⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ Φ Φ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ Ψ Ψ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦

Solve each time step using SuperLU For linear problem,only need to form LU decomposition

• nce and do a

back-substitution each time step. Note that stream function and vorticity are solved together

Tilting of a Plasma Column

1 1

[2/ ( )] ( )cos , 1 ( 1/ )cos , 1 ( ) kJ k J kr r r r r J k θ ψ θ < ⎧ = ⎨ − > ⎩ =

Initial Condition: Give small perturbation and evolve in time

Stream function and vorticity at final time Flux (top) and current (bottom) at initial and final times

Tilting of a Plasma Column-cont

Converged (in time) growth rate the same for N=30,40 out to 6 decimal places

Boundary Conditions

x y

x y xx xy yy

φ φ φ φ φ φ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ Φ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ If boundary aligned with x axis: Homogeneous Dirichlet Homogeneous Neumann

x xx

φ φ φ = = =

y xy

φ φ = =

Non-Rectangular Boundary Conditions

At each boundary point, define a local normal and curvature: ˆ ˆ ˆ / n n dt dl κ ≡ ⋅ We then define a new set of basis functions for that boundary node that are locally aligned with the boundary and curvature Old basis functions Ignoring curvature with curvature

2 2 2 2 2 2 2 2 2 2

1 2 ( ) 2 2 2

x y y x y x x n y x x y x y x y y t x x y y xx nn x y x y x y xy nt y x y x yy tt

n n n n n n n n n n n n n n n n n n n n n n n n n n n n ν μ κ κ κ ν μ κ κ κ ν μ ν μ ν μ ν μ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − − − − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ i

34

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

35

2-Fluid MHD Equations:

( )

2 2 2 2 2

( ) ( ) 3 3 2 2 3 3 5 2 2 2 1 3 2

i e e e e GV e H e e i i i i

n n t t nM p t p p p J Q t p p h ne p p n ne p p V Q t n λ μ η η μ

Δ Δ

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

• ∇

+ ∇ = × − + ∇ ∂ + × = + ∂ ⎛ ⎞ + ∇ × −∇ − ∇ ⎡ ⎤ ∇ − ∇ ⎢ ⎥ ⎣ ⎦ ∇ = − ∇ + + − ∇ + ⎜ ⎟ ∂ ⎝ ⎠ ∂ ⎛ ⎞ + ∇ = − ∇ + ∇ −∇ − ⎜ ⎟ ∂ ⎝ ⎠ V B E B A J B V V V J B V E V Π J B J B J V V q V q J V i i i i i i i i

2 R

• es

fl isti uid ve MHD terms

2 2 2

R R U R χ ϕ ϕ ω

− ⊥

= ∇ ×∇ + ∇ + ∇ V

2

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

8 scalar 3D variables: , , , , , , ,

e i

f U n p p ψ χ ω

No further approximations! Solve these as faithfully as possible in realistic geometry ϕ Z R

36

t = 1 t = 16 t = 24 t = 32 t = 40 t = 8

“GEM” reconnection test problem1 for 2-fluid MHD

• Starts like resistive MHD
• Dramatic change in configuration for t > 20

In-plane current density contours at different times

• 1J. Birn, et al, J. Geophys.
• Res. 106 (2001) 3715

37

2-fluid reconnection requires high resolution for convergent results

• Note sudden transition where

velocity abruptly increases

• These calculations used a

hyperviscosity term in Ohm’s law proportional to h2 . . . required for a stable calculation

• Reconnected Flux at t=40

converges as 1/h5

38

vorticity field velocity divergence

• ut-of-plane current
• ut-of-plane magnetic field

2F GEM Reconnection snapshot at time of maximum velocity t ~ 32 (2002 nodes) Energy conservation to 1 in 103, flux conservation exact, symmetry preserved even though triangles were not arranged symmetrically. midplane

39

Midplane Current density collapses to the width

• f 1-3 triangular elements

t=32 time of previous contour plot ( note sudden collapse at t=23+)

40

Midplane electric field before and after transition

t=20

X-Midplane

10 15 20 25 30 35

• 0.01

0.00 0.01 0.02 0.03 0.04

t=30

X-Midplane

10 15 20 25 30 35

• 0.2

0.0 0.2 0.4

Total Electric Field Resistivity J x B V x B Hyper-resistivity

( )

2 2

1 ˆ

e H

J B p h J V E ne B J z η λ × − ⎡ ⎤ = − ∇ + + ⎢ ⎥ ⎣ ⎦ × −∇

• i
• Reconnection rate:

41

Blowup showing electric field after transition

X-Midplane

20 21 22 23 24 25

• 0.2

0.0 0.2 0.4

Total Electric Field Resistivity J x B V x B Hyper-resistivity

( )

2 2

1 ˆ

e H

J B p h J V E ne B J z η λ × − ⎡ ⎤ = − ∇ + + ⎢ ⎥ ⎣ ⎦ × −∇

• i
• ne triangle

Hyper-resistivity coefficient must be large enough that current density collapse is limited to 1-2 triangles: reason for factor h2

t=30

X-Midplane

10 15 20 25 30 35 0.2 0.0 0.2 0.4

42

The split-implicit time advance

43

2 2 2 2 2

u v c u u t x c v u t x c t x ∂ ∂ ⎫ = ⎪ ∂ ∂ ⎪ ∂ ∂ − = ⎬ ∂ ∂ ∂ ∂ ⎪ = ⎪ ∂ ∂ ⎭

Consider a simple 1-D Hyperbolic System of Equations (Wave Equation)

44

2 2 2 2 2

u v c u u t x c v u t x c t x ∂ ∂ ⎫ = ⎪ ∂ ∂ ⎪ ∂ ∂ − = ⎬ ∂ ∂ ∂ ∂ ⎪ = ⎪ ∂ ∂ ⎭

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 θ θ δ δ δ θ θ δ δ δ

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

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

Consider a simple 1-D Hyperbolic System of Equations (Wave Equation) Implicit Centered Difference: Stable for θ ≥ ½ 2nd order for θ = ½

45

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

46

2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S S ⎡ ⎤ + − − ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ − + − − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − + − − ⎢ ⎥ ⎢ ⎥ − + − − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − + − − ⎢ ⎥ ⎢ ⎥ − + − ⎢ ⎥ ⎣ ⎦ ⎢ ⎥ → ⎢ ⎥ − ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦

2N N N Substitution takes us from having to invert a 2N x 2N system that has large off-diagonal elements to sequentially inverting a N x N system that is diagonally dominant + the identity matrix. Mathematically equivalent same answers!

Matrices to be inverted

1 c t S x θ δ δ =

• Coupled system

Un-coupled system

47

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

1 1 (1 ) (1 )

n n n x j x j x j n n n n j j x j x j

S u S u S v v v S u u θ δ θ θ δ δ θδ θ δ

+ + + + + + +

⎡ ⎤ ⎡ ⎤ − = + − + ⎣ ⎦ ⎣ ⎦ ⎡ ⎤ = + + − ⎣ ⎦

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

2

n n n n x j j j j n n n x j j j

u u u u v v v tc S 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 δ θ θ θ δ δ δ δ δ θ θ δ

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

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

Rewrite using standard finite-difference notation: Operator to invert

48

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

49

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

50

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:

51

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

Summary

• Finite elements with C1 continuity can be used with differential equations up

to 4th order.

• Hermite cubic element is an attractive C1 element in 1D

– 2 DOF at each node

• Q18 is an attractive C1 element in 2D

– All DOF defined at nodes very compact compared to other elements – Shown to be very accurate – Curved domains can be accommodated with normal vector and curvature – High accuracy has been verified in demanding tests

• 3D Element can be made as the tensor product of Q18 and Hermite cubic
• Fits naturally into implicit method for MHD in tokamak (next lecture)

53