A globally divergence-free discontinuous galerkin method for - PowerPoint PPT Presentation

A globally divergence-free discontinuous galerkin method for induction and related equations Praveen Chandrashekar praveen@math.tifrbng.res.in Center for Applicable Mathematics Tata Institute of Fundamental Research Bangalore-560065, India http://cpraveen.github.io Int. Conf. on Current Trends in Theoretical and Computational Differential Equations with Applications, SAU, New Delhi 1–5 December 2017 Supported by Airbus Foundation Chair at TIFR-CAM, Bangalore http://math.tifrbng.res.in/airbus-chair 1 / 54

Joint work with Dinshaw Balsara Univ. of Notre Dame 2 / 54

Maxwell Equations Linear hyperbolic system ∂ B ∂ D ∂t + ∇ × E = 0 , ∂t − ∇ × H = − J B = magnetic flux density D = electric flux density E = electric field H = magnetic field J = electric current density µ, ε ∈ R 3 × 3 symmetric B = µ H , D = ε E , J = σ E ε = permittivity tensor µ = magnetic permeability tensor σ = conductivity ∂ρ ∇ · B = 0 , ∇ · D = ρ (electric charge density) , ∂t + ∇ · J = 0 3 / 54

Ideal MHD equations Nonlinear hyperbolic system Compressible Euler equations with Lorentz force ∂ρ ∂t + ∇ · ( ρ v ) = 0 ∂ρ v ∂t + ∇ · ( pI + ρ v ⊗ v − B ⊗ B ) = 0 ∂E ∂t + ∇ · (( E + p ) v + ( v · B ) B ) = 0 ∂ B ∂t − ∇ × ( v × B ) = 0 Magnetic monopoles do not exist: = ⇒ ∇ · B = 0 4 / 54

Model problem ∂ B ∂t + ∇ × E = − M Divergence evolves according to ∂ ∂t (div B ) + ∇ · M = 0 (1) In 2-D ∂B x ∂t + ∂E z ∂B y ∂t − ∂E z ∂y = − M x , ∂x = − M y 5 / 54

Model problem In MHD, B represents the magnetic field and M = 0 ∂ B ∂t + ∇ × E = 0 , E = − v × B Magnetic monopoles do not exist: = ⇒ ∇ · B = 0 If ∇ · B = 0 at t = 0 then ∂ ∂t ( ∇ · B ) + ∇ · ∇ × E = 0 = ⇒ ∇ · B = 0 for t > 0 In 2-D, the induction equation can be written as ∂B x ∂t + ∂E ∂B y ∂t − ∂E ∂y = 0 , ∂x = 0 , E = v y B x − v x B y 6 / 54

Some existing methods for MHD Exactly divergence-free methods • Constrained transport ([1] Evans & Hawley (1989)) ◮ ∇ · B = 0 implies B = ∇ × A ◮ Evolve A forward in time ◮ Compute B from A • Divergence-free reconstruction ([2] Balsara (2001)) • Globally divergence-free scheme ([3] Li et al. (2011)) Approximate methods • Using Godunov’s symmetrized version of MHD [4] (Powell [5], CK [6]) • Divergence cleaning methods (Dedner et al. [7]) 7 / 54

Approximation of magnetic field When dealing with problems where the vector field B must be divergence-free, it is natural to look for solutions in H (div , Ω) which is defined as H (div , Ω) = { B ∈ L 2 (Ω) : div( B ) ∈ L 2 (Ω) } To approximate functions in H (div , Ω) on a mesh T h , we need the following compatibility condition. Theorem (See [8], Proposition 3.2.2) Let B h : Ω → R d be such that 1 B h | K ∈ H 1 ( K ) for all K ∈ T h 2 for each common face F = K 1 ∩ K 2 , K 1 , K 2 ∈ T h , the trace of normal component n · B h | K 1 and n · B h | K 2 is the same. Then B h ∈ H ( div , Ω) . Conversely, if B h ∈ H ( div , Ω) and (1) holds, then (2) is also satisfied. 8 / 54

Approximation of magnetic field P k ( x ) , P k ( y ) : 1-D polynomials of degree at most k wrt the variables x , y respectively. Q r,s ( x, y ) : tensor product polynomials of degree r in the variable x and degree s in the variable y , i.e., Q r,s ( x, y ) = span { x i y j , 0 ≤ i ≤ r, 0 ≤ j ≤ s } For k ≥ 0 , the Raviart-Thomas space of vector functions is defined as RT k = Q k +1 ,k × Q k,k +1 , dim( RT k ) = 2( k + 1)( k + 2) Two consequences: • For any B h ∈ RT k , we have div( B h ) ∈ Q k,k ( x, y ) =: Q k ( x, y ) 9 / 54

Approximation of magnetic field • The restriction of B h = ( B h x , B h y ) to a face is a polynomial of degree k , i.e., B h B h x ( ± ∆ x/ 2 , y ) ∈ P k ( y ) , y ( x, ± ∆ y/ 2) ∈ P k ( x ) For doing the numerical computations, it is useful to map each cell to a reference cell. { ξ i , 0 ≤ i ≤ k + 1 } = Gauss-Lobatto-Legendre (GLL) nodes { ˆ ξ i , 0 ≤ i ≤ k } = Gauss-Legendre (GL) nodes Let φ i and ˆ φ i be the corresponding 1-D Lagrange polynomials. Then the magnetic field is given by k +1 k k k +1 � � ( B x ) ij φ i ( ξ )ˆ � � ( B y ) ij ˆ B h B h x ( ξ, η ) = φ j ( η ) , y ( ξ, η ) = φ i ( ξ ) φ j ( η ) i =0 j =0 i =0 j =0 10 / 54

Approximation of magnetic field B y B x Location of dofs of Raviart-Thomas polynomial for k = 0 B x B y Location of dofs of Raviart-Thomas polynomial for k = 1 11 / 54

Approximation of magnetic field B y B x Location of dofs of Raviart-Thomas polynomial for k = 2 Our choice of nodes ensures that the normal component of the magnetic field is continuous on the cell faces. We have the error estimates on Cartesian meshes [9], [10] � B − B h � L 2 (Ω) ≤ Ch k +1 | B | H k +1 (Ω) � div( B ) − div( B h ) � L 2 (Ω) ≤ Ch k +1 | div( B ) | H k +1 (Ω) 12 / 54

Construction of B h from moments e + 2 y 3 The edge moments are given by � B h x φ d y ∀ φ ∈ P k ( y ) e ∓ x e − e + C x x and � B h y φ d x ∀ φ ∈ P k ( x ) e ∓ y 0 1 e − y The cell moments are given by � B h x ψ d x d y ∀ ψ ∈ ∂ x Q k ( x, y ) := Q k − 1 ,k ( x, y ) C 13 / 54

Construction of B h from moments and � B h y ψ d x d y ∀ ψ ∈ ∂ y Q k ( x, y ) := Q k,k − 1 ( x, y ) C Note that dim P k ( x ) = dim P k ( y ) = k + 1 and dim ∂ x Q k ( x, y ) = dim ∂ y Q k ( x, y ) = k ( k + 1) so that we have in total 4( k + 1) + 2 k ( k + 1) = 2( k + 1)( k + 2) = dim RT k The moments on the edges e ∓ x uniquely determine the restriction of B h x on those edges, and similarly the moments on e ∓ y uniquely determine the restriction of B h y on the corresponding edges. This ensures continuity of the normal component of B h on all the edges. 14 / 54

Theorem If all the moments are zero for any cell C , then B h ≡ 0 inside that cell. Proof: The edge moments being zero implies that e ∓ e ∓ B h B h x ≡ 0 on and y ≡ 0 on x y Now take ψ = ∂ x φ for some φ ∈ Q k in the cell moment equation of B h x and perform an integration by parts ∂B h � � � x B h B h − ∂x φ d x d y − x φ d y + x φ d y = 0 e + e − C x x and hence ∂B h � x ∂x φ d x d y = 0 ∀ φ ∈ Q k C Since ∂B h ∂x ∈ Q k , this implies that ∂B h ∂x ≡ 0 and hence B h x ≡ 0 . Similarly, x x we conclude that B h y ≡ 0 . 15 / 54

Theorem Let B h ∈ RT k satisfy the moments � � B h x φ d y = B x φ d y ∀ φ ∈ P k ( y ) (2) e ∓ e ∓ x x � � B h y φ d x = B y φ d x ∀ φ ∈ P k ( x ) (3) e ∓ e ∓ y y � � B h x ψ d x d y = B x ψ d x d y ∀ ψ ∈ ∂ x Q k ( x, y ) (4) C C � � B h y ψ d x d y = B y ψ d x d y ∀ ψ ∈ ∂ y Q k ( x, y ) (5) C C for a given vector field B ∈ H ( div , Ω) . If div ( B ) ≡ 0 then div ( B h ) ≡ 0 . 16 / 54

Proof: We choose ψ = ∂ x φ and ψ = ∂ y φ for some φ ∈ Q k ( x, y ) respectively in the two cell moment equations (4), (5). Adding these two equations together, we get � � ( B h x ∂ x φ + B h y ∂ x φ ) d x d y = ( B x ∂ x φ + B y ∂ y φ ) d x d y C C Performing integration by parts on both sides � � � � div( B h ) φ d x d y + φ B h · n d s = − − div( B ) φ d x d y + φ B · n d s C ∂C C ∂C Note that φ restricted to ∂C is a one dimensional polynomial of degree k and the edge moments of B h and B agree with one another by equations (2), (3). Hence we get � � div( B h ) φ d x d y = div( B ) φ d x d y ∀ φ ∈ Q k ( x, y ) C C 17 / 54

If div( B ) ≡ 0 then � div( B h ) φ d x d y = 0 ∀ φ ∈ Q k ( x, y ) C Since div( B h ) ∈ Q k ( x, y ) this implies that div( B h ) ≡ 0 everywhere inside the cell C . Remark : The proof makes use of integration by parts for which the quadrature must be exact. The integrals involving B h can be evaluated exactly using Gauss quadrature of sufficient accuracy. This is not the case for the integrals involving B since it can be an arbitrary nonlinear function. When div( B ) = 0 , we have B = ( ∂ y Φ , − ∂ x Φ) for some smooth function Φ . We can approximate Φ by Φ h ∈ Q k +1 and compute the projections using ( ∂ y Φ h , − ∂ x Φ h ) in which case the integrations can be performed exactly. 18 / 54

Example: RT 0 B h B h x ( x, y ) = a 0 + a 1 x, y ( x, y ) = b 0 + b 1 y In this case we have only the edge moments. The polynomial test function spaces needed to specify the edge moments are P 0 ( x ) = span { 1 } , P 0 ( y ) = span { 1 } and the four moments corresponding to the four faces are 1 1 � � 2 2 B h B h x ( − 1 / 2 , y ) d y = α 1 x (1 / 2 , y ) d y = α 2 − 1 − 1 2 2 1 1 � � 2 2 B h B h y ( x, − 1 / 2) d x = β 1 y ( x, 1 / 2) d x = β 2 − 1 − 1 2 2 The solution is given by a 0 = 1 2 ( α 1 + α 2 ) , a 1 = α 2 − α 1 b 0 = 1 2 ( β 1 + β 2 ) , b 1 = β 2 − β 1 19 / 54

Example: RT 1 a 0 + a 1 x + a 2 y + a 3 xy + a 4 ( x 2 − 1 12 ) + a 5 ( x 2 − 1 B h x ( x, y ) = 12 ) y b 0 + b 1 x + b 2 y + b 3 xy + b 4 ( y 2 − 1 12 ) + b 5 x ( y 2 − 1 B h y ( x, y ) = 12 ) The polynomial test function spaces needed to specify the moments (2)-(5) are P 1 ( x ) = span { 1 , x } , P 1 ( y ) = span { 1 , y } ∂ x Q 1 ( x, y ) = span { 1 , y } , ∂ y Q 1 ( x, y ) = span { 1 , x } 20 / 54