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 Today Tomorrow • Galerkin method • Split implicit time C 1 finite elements • differencing for MHD • Examples • Accuracy, spectral pollution, and • Split implicit time representation of the differencing 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]: ⎡ ⎤ ⎛ ⎞ d d − + = ⎜ ( ) ⎟ ( ) ( ) ( ) (1) p x q x U x f x ⎢ ⎥ ⎝ ⎠ ⎣ ⎦ dx dx We choose a function space H 1 , and multiply (1) by an arbitrary member of that function space, and integrate over the domain: ⎡ ⎤ b ⎛ ⎞ b d d ∫ ∫ − + = ( ) ⎜ ( ) ⎟ ( ) ( ) ( ) ( ) (2 ) V x p x q x U x dx V x f x a ⎢ ⎥ ⎝ ⎠ ⎣ ⎦ dx dx a a Or, after integrating by parts (assuming the boundary term vanishes): ⎡ ⎤ b ⎛ ⎞ + b dV dU ∫ ∫ = ⎜ ⎟ ( ) ( ) ( ) ( ) (2 ) p x q x VU dx V x f x b ⎢ ⎥ ⎝ ⎠ ⎣ ⎦ dx dx a a If we require that Eq. (2) hold for every function V(x) in the function space H 1 , it is called the weak form of Eq. (1). Note that only first derivatives occur in (2b), whereas 2 nd derivatives occur in (1)
Galerkin Finite Element Method The Galerkin finite element method is a discretization of the weak form of the equation. To apply, we chose a finite dimensional subspace S of the infinite dimensional Hilbert space H 1 , and require that Eq. (2) hold for every function in S. N = ∑ ϕ known basis functions that span S ( ) ( ) U x a x j j a j are amplitudes to be solved for = j 1 = ϕ Letting Equation (2b) becomes ( ) ( ) V x x i ϕ ⎡ ⎤ ⎛ ϕ ⎞ + b b N d ∑ ∫ d ∫ ϕ ϕ = ϕ = j ⎢ ( ) i ( ) ⎥ ( ) ( ) 1, a ⎜ p x ⎟ q x dx x f x i N j i j i ⎝ ⎠ ⎣ ⎦ dx dx = j 1 a a Or, = M a b ij j i ϕ ⎡ ⎤ ⎛ ϕ ⎞ b d d ∫ = + ϕ ϕ j ⎢ ( ) i ( ) ⎥ M ⎜ p x ⎟ q x dx ij i j ⎝ ⎠ ⎣ ⎦ dx dx a b ∫ = ϕ ( ) ( ) b x f x dx i i a In the finite element method, the basis functions are low order polynomials
Simple Example: Linear Elements in 1D ϕ ⎡ ⎤ ⎛ ϕ ⎞ + b b = ∑ N d N ∑ ∫ d ∫ ϕ ϕ = ϕ = ϕ j ⎢ ( ) i ( ) ⎥ ( ) ( ) 1, ( ) ( ) a ⎜ p x ⎟ q x dx x f x i N U x a x j i j i j j ⎝ ⎠ ⎣ ⎦ dx dx = = j 1 j 1 a a = = = , (constants) ( ) p p q q f f x 0 0 Linear ϕ = δ = , at x x ϕ ϕ + i i j j ϕ − ( ) 1 ( ) elements only i x x 1 ( ) x i i overlap with linear: x i = ih their nearest 1 h=L/N neighbors = ∫∫ ϕ b fdx 0 x i-1 x i x i+1 L i i − ⎡ ⎤ ⎡ ⎤ 2 1 0 0 4 1 0 0 ⎢ ⎥ ⎢ ⎥ − − 1 2 1 0 1 4 1 0 1 ⎢ ⎥ h ⎢ ⎥ ∫∫ ∫∫ ∇ ϕ ∇ ϕ = = = = i p dx K p h q v v dx A q ⎢ ⎥ − − ⎢ ⎥ 0 , 0 i j i j 0 i j i j , 0 0 1 2 1 6 0 1 4 1 ⎢ ⎥ ⎢ ⎥ − ⎣ 0 0 1 2 ⎦ ⎣ 0 0 1 4 ⎦ “stiffness matrix” “mass matrix” Leads to sparse, diagonally = M a b dominant matrices: M i,j = K i,j + A i,j ij j i
Note: We were able to solve a 2 nd order ODE with linear elements! ⎡ ⎤ ⎛ ⎞ b d d − + = ∫ ϕ ⎜ ( ) ⎟ ( ) ( ) ( ) p x q x U x f x = ∑ ⎢ ⎥ N ( ) dx x ϕ ⎝ ⎠ ⎣ ⎦ i dx dx ( ) ( ) U x a x j j a = j 1 ϕ ⎡ ⎤ ⎛ ϕ ⎞ + b b N d ∑ ∫ d ∫ ϕ ϕ = ϕ = j ⎢ ( ) i ( ) ⎥ ( ) ( ) 1, a ⎜ p x ⎟ q x dx x f x i N j i j i ⎝ ⎠ ⎣ ⎦ dx dx = j 1 a a ϕ ϕ + ϕ − ( ) 1 ( ) i x x 1 ( ) x i i 1 0 x i-1 x i x i+1 L = M a b ij j i
Extension to Higher Order Equations Suppose we have a 4 th order equation: 4 = ∑ N d U x = ϕ ( ) ( ) ( ) ( ) f x U x a x 4 j j dx = j 1 Form the weak form: b b 4 d ∫ ∫ ϕ = ϕ ( ) ( ) ( ) ( ) x U x dx x f x dx i i 4 dx a a ϕ b b 3 d d U dx ∫ ∫ − = ϕ i ( ) ( ) x f x dx Integrate by parts twice: (assumes the 3 i dx dx a a boundary ϕ b b 2 2 d d U dx ∫ ∫ terms vanish) = ϕ i ( ) ( ) x f x dx 2 2 i dx dx a a Application of the Galerkin method to a 4 th order equation requires that the second derivative of the basis functions be defined everywhere. For a 2 nd order equation, only the first derivative need be defined.
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 H 2 ϕ ϕ 2 ( ) 1 ( ) j x j x The Hermite Cubic elements associate two basis functions with each node j, defined on the interval x j-1 < x < x j+1 as: ( ) ( ) 2 ϕ = − + ( ) 1 2 1 y y y j-1 j j+1 j-1 j j+1 1 ( ) 2 ϕ = − ( ) 1 y y y These have the property that: 2 − ⎛ ⎞ 1. In the interval x j < x < x j+1 , a complete x x ϕ = ϕ j ( ) x ⎜ ⎟ cubic can be represented by: 1, j 1 ⎝ ⎠ h ϕ ϕ ϕ ϕ ( ), ( ), ( ), ( ), x x x x + + − 1, j 2, j 1, j 1 2, j 1 ⎛ ⎞ x x ϕ = ϕ j ( ) x h ⎜ ⎟ 2. The first derivatives vanish at the 2, 2 j ⎝ ⎠ h element boundaries 3. The second derivatives are defined everywhere
Hermite Cubic Finite Elements ϕ ϕ 2 ( ) 1 ( ) j x j x N ∑ = ⎡ ϕ + ϕ ⎤ ( ) ( ) ( ) U x a x b x ⎣ ⎦ 1, 2, j j j j = 1 j In the interval [ j , j+1 ] , the function is written as: = ϕ + ϕ + ϕ + ϕ ( ) ( ) ( ) ( ) ( ) U x a x b x a x b x j-1 j j+1 j-1 j j+1 + + + + 1, 2, 1 1, 1 1 2, 1 j j j j j j j j ( ) ( ) 2 = + + − − + − 3 2 3 a b x a hb a hb x h Simple physical interpretation: + + 1 1 j j j j j j ( ) ( ) + + − + 3 2 2 a hb a hb x h + + a j is value of function at node j j j j 1 j 1 = + + + b j is value of derivative at node j 2 3 c c x c x c x 0 1 2 3 ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ 1 0 0 0 a c Complete cubic 0 j ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ polynomial determined 0 1 0 0 b c ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = 1 j i by values of function ⎢ ⎥ ⎢ ⎥ ⎢ − − − ⎥ 2 2 3 / 2 / 3 / 1/ a c h h h h and derivative at + 1 2 j ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − endpoints 3 2 3 2 ⎣ ⎦ ⎣ ⎦ b 2 / 1/ 2 / 1/ ⎣ ⎦ c h h h h + 1 3 j
Hermite Cubic Finite Elements-2 N ∑ = ⎡ ϕ + ϕ ⎤ ϕ ( ) ( ) ( ) 1 ( ) U x a x b x ⎣ ⎦ j x 1, 2, j j j j = 1 j In the interval [ j , j+1 ] , the function is written as: ϕ 2 ( ) j x = ϕ + ϕ + ϕ + ϕ ( ) ( ) ( ) ( ) ( ) U x a x b x a x b x + + + + 1, 2, 1 1, 1 1 2, 1 j j j j j j j j ( ) ( ) 2 = + + − − + − 3 2 3 a b x a hb a hb x h + + 1 1 j j j j j j ϕ 1 ( ) ( ) ( x ) + + − + 3 + 2 2 1 j a hb a hb x h + + j j j 1 j 1 = + + + 2 3 c c x c x c x 0 1 2 3 ϕ 1 ( ) x Complete cubic in each interval! + 2 j j j+1
Comparison of elements and derivatives 1 st Hermite 2 nd Hermite Linear cubic cubic elements ϕ ( ) x ϕ ′ ( ) x j-1 j j+1 ϕ ′′ ( ) x not defined j-1 j j+1 j-1 j j+1
Types of Finite Elements in 2D Shape Order Continuity = ∑ φ C 0 : continuous i k ( , ) x y C x y ik across element + < i k M ≥ boundaries , 0 rectangle i k C 1 : continuous Order M if complete first derivatives polynomial to that across element order boundaries triangle DG: discontinuous Galerkin
Element Shape Rectangle: Triangle: • natural for structured mesh • natural for unstructured mesh • allows tensor product forms • easier to fit complex shapes = ∑ φ i k ( , ) x y a x y • can easily refine locally ik ≤ ≤ 0 , i k t bi-linear: t=1, bi-quadratic: t=2, etc •cannot refine locally without hanging nodes
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 h M+1 This follows directly from a local Taylor series expansion: ⎡ ∂ φ ⎤ k M k 1 ∑∑ − + φ = − − + 1 l k l M ( , ) ( ) ( ) ( ) x y ⎢ ⎥ x x y y O h − ∂ ∂ − 0 0 l k l !( )! ⎣ ⎦ l k l x y = = 0 0 k l x , y 0 0 Thus, linear elements will be O( h 2 ) quadratic elements will be O( h 3 ) cubic elements will be O( h 4 ) quartic elements will be O( h 5 ) complete quintic elements will be O( h 6 )
Recommend
More recommend
