Divergence-free discontinuous Galerkin method for MHD and Maxwell’s equations Praveen Chandrashekar praveen@math.tifrbng.res.in Center for Applicable Mathematics Tata Institute of Fundamental Research Bangalore-560065, India http://cpraveen.github.io In collaboration with Dinshaw Balsara, Arijit Hazra, Rakesh Kumar Zurich Colloquium in Applied Mathematics 10 April 2019 1 / 38
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 2 / 38
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 3 / 38
Divergence constraint Rotated shock tube Powell 0 . 1 ∂ B ∂t + ∇ × E = 0 0 . 0 − 0 . 1 Relative error on B � ∇ · ∂ B cleaning 0 . 1 ∂t + ∇ · ∇× E = 0 � �� � 0 . 0 =0 − 0 . 1 ∂ ∂t ∇ · B = 0 Athena 0 . 1 0 . 0 If − 0 . 1 ∇ · B = 0 at t = 0 − 0 . 4 − 0 . 2 0 . 0 0 . 2 0 . 4 x then Guillet et al., MNRAS 2019 ∇ · B = 0 for t > 0 Discrete div-free = ⇒ positivity (Kailiang) 4 / 38
Objectives • Divergence-free schemes for Maxwell’s and compressible MHD ◮ Cartesian grids at present ◮ Divergence-free reconstruction (BDM) ◮ Divergence preserving schemes (RT) • High order accurate ◮ discontinuous-Galerkin • Upwind-type schemes based on Riemann solvers • Non-oscillatory schemes for MHD • Explicit time stepping 5 / 38
Some existing methods Exactly divergence-free methods • Yee scheme (Yee (1966)) • Projection methods (Brackbill & Barnes (1980)) • Constrained transport (Evans & Hawley (1989)) • Divergence-free reconstruction (Balsara (2001)) • Globally divergence-free scheme (Li et al. (2011), Fu et al, (2018)) Approximate methods • Locally divergence-free schemes (Cockburn, Li & Shu (2005)) • Godunov’s symmetrized version of MHD (Powell, Gassner et al., C/K) • Divergence cleaning methods (Dedner et al.) 6 / 38
Approximation of magnetic field If ∇ · B = 0 , it is natural to look for approximations in H ( div, Ω) = { B ∈ L 2 (Ω) : div( B ) ∈ L 2 (Ω) } To approximate functions in H ( div, Ω) on a mesh T h with piecewise polynomials, we need B · n continuous across element faces. 7 / 38
DG scheme based on divergence-free reconstruction of 2-D vector fields 8 / 38
Approximation spaces B · n on the faces must be continuous: approximate B · n by polynomial functions P k on the faces. 2 3 [ ∆ x , ∆ y ] → [ − 1 2 , + 1 2 ] × [ − 1 2 , + 1 2 ] B + y ( ξ ) ξ = x − x c η = y − y c ∆ x , ∆ y k B − B + x ( η ) x ( η ) C � B ± a ± x ( η ) = j φ j ( η ) j =0 k � B ± b ± y ( ξ ) = j φ j ( ξ ) B − y ( ξ ) j =0 0 1 9 / 38
Approximation spaces φ j : [ − 1 2 , 1 2 ] → R are mutually orthogonal polynomials; hence 1 1 � � 2 2 B ± x ( η ) d η = a ± B ± y ( ξ ) d ξ = b ± 0 , 0 − 1 − 1 2 2 Given the data B ± x , B ± y on the faces, we want to recontruct a divergence-free vector field ( B x ( ξ, η ) , B y ( ξ, η )) inside the cell. � � Necessary condition: B · n = ∇ · B = 0 ∂C C Writing this condition in the reference coordinates yields 1 1 � � 2 2 ( B + x ( η ) − B − ( B + y ( ξ ) − B − ∆ y x ( η )) d η + ∆ x y ( ξ )) d ξ = 0 − 1 − 1 2 2 or � B · n = ( a + 0 )∆ y + ( b + 0 − a − 0 − b − 0 )∆ x = 0 ∂C 10 / 38
Raviart-Thomas (RT) polynomials Space of tensor product polynomials Q r,s = span { ξ i η j : 0 ≤ i ≤ r, 0 ≤ j ≤ s } For k ≥ 0 , approximate the vector field inside the cells by tensor product polynomials k +1 k � � B x ( ξ, η ) = a ij φ i ( ξ ) φ j ( η ) ∈ Q k +1 ,k i =0 j =0 k k +1 � � B y ( ξ, η ) = b ij φ i ( ξ ) φ j ( η ) ∈ Q k,k +1 i =0 j =0 RT( k ) = Q k +1 ,k × Q k,k +1 , dim RT( k ) = 2( k + 1)( k + 2) div RT( k ) ∈ Q k,k 11 / 38
Brezzi-Douglas-Marini (BDM) polynomials The BDM polynomials are of the form (Brezzi & Fortin, Section III.3.2) BDM( k ) = ( P k ) 2 + r ∇ × ( x k +1 y ) + s ∇ × ( xy k +1 ) dim BDM( k ) = ( k + 1)( k + 2) + 2 , div BDM( k ) ∈ P k − 1 In reference coordinates, we take the polynomials to be of the form k a ij φ i ( ξ ) φ j ( η ) + r ( ξ k +1 − M k +1 ) + s ( k + 1) ξη k � B x ( ξ, η ) = ∆ y ∆ y i,j =0 i + j ≤ k k − s ( η k +1 − M k +1 ) b ij φ i ( ξ ) φ j ( η ) − r ( k + 1) ξ k η � B y ( ξ, η ) = ∆ x ∆ x i,j =0 i + j ≤ k 12 / 38
Brezzi-Douglas-Fortin-Marini (BDFM) polynomials The BDFM polynomials are of the form BDFM( k ) = ( P k +1 ) 2 \ (0 , x k +1 ) \ ( y k +1 , 0) dim BDFM( k ) = ( k + 2)( k + 3) − 2 , div BDFM( k ) ∈ P k In reference cell coordinates, the polynomials can be written as k +1 k − i k +1 k − j � � � � B x ( ξ, η ) = a ij φ i ( ξ ) φ j ( η ) , B y ( ξ, η ) = b ij φ i ( ξ ) φ j ( η ) i =0 j =0 j =0 i =0 Remark: Raviart-Thomas space is the largest space considered here, and in fact we have BDM( k ) ⊂ BDFM( k ) ⊂ RT( k ) 13 / 38
Conditions for the reconstruction For all the three polynomial spaces B x ( ± 1 B y ( ξ, ± 1 2 , η ) ∈ P k and 2 ) ∈ P k Match B ( ξ, η ) to B · n on the four faces B x ( ± 1 2 , η ) = B ± ∀ η ∈ [ − 1 2 , + 1 B y ( ξ, ± 1 2 ) = B ± ∀ ξ ∈ [ − 1 2 , + 1 x ( η ) 2 ] , y ( ξ ) 2 ] In addition, we have to make the vector field divergence-free 1 ∂B x ∂ξ + 1 ∂B y ∀ ξ, η ∈ [ − 1 2 , + 1 ∇ · B = ∂η = 0 , 2 ] ∆ x ∆ y Questions : • Do we have enough equations ? (Not always) • Can we solve them ? (Yes) Approach : • All three spaces contain P k , necessary to obtain h k +1 ’th order accuracy. • Throw away any basis function outside P k without affecting the order of accuracy. 14 / 38
Summary of reconstruction • All three spaces lead to same solution: BDM (k) • k = 0 , 1 , 2 : reconstruction determined by face solution alone • k ≥ 3 : Require additional information to solve reconstruction problem 15 / 38
Reconstruction at degree k = 3 The face solution is a polynomial of degree three and is of the form Right/left face B ± x ( η ) = a ± 0 φ 0 ( η ) + a ± 1 φ 1 ( η ) + a ± 2 φ 2 ( η ) + a ± 3 φ 3 ( η ) Top/bottom face y ( η ) = b ± 0 φ 0 ( ξ ) + b ± 1 φ 1 ( ξ ) + b ± 2 φ 2 ( ξ ) + b ± B ± 3 φ 3 ( ξ ) 16 / 38
Reconstruction using BDFM(3) The vector field has the form B x ( ξ, η ) = a 00 + a 10 ξ + a 01 η + a 20 ( ξ 2 − 12 ) + a 11 ξη + a 02 ( η 2 − 1 1 12 )+ a 30 ( ξ 3 − 20 ξ ) + a 21 ( ξ 2 − 12 ) η + a 12 ξ ( η 2 − 12 ) + a 03 ( η 3 − 3 1 1 3 20 η )+ a 40 ( ξ 4 − 14 ξ 2 + 560 ) + a 13 ξ ( η 3 − 20 η ) + a 31 ( ξ 3 − 3 3 3 3 20 ξ ) η + a 22 ( ξ 2 − 12 )( η 2 − 1 12 ) ∈ P 4 \ { y 4 } 1 B y ( ξ, η ) = b 00 + b 10 ξ + b 01 η + b 20 ( ξ 2 − 12 ) + b 11 ξη + b 02 ( η 2 − 1 1 12 )+ b 30 ( ξ 3 − 20 ξ ) + b 21 ( ξ 2 − 12 ) η + b 12 ξ ( η 2 − 12 ) + b 03 ( η 3 − 3 1 1 3 20 η )+ b 22 ( ξ 2 − 12 )( η 2 − 12 ) + b 13 ξ ( η 3 − 20 η ) + b 31 ( ξ 3 − 1 1 3 3 20 ξ ) η + b 04 ( η 4 − 14 η 2 + 560 ) ∈ P 4 \ { x 4 } 3 3 which has a total of 28 coefficients . The divergence-free condition gives 9 equations and matching the face solution gives 16 equations for a total of 26 equations . 17 / 38
Reconstruction using BDFM(3) We can first solve for some of the coefficients completely in terms of the face solution 1 1 ∆ x 0 + a + ( b + 1 1 ∆ y ( a − 1 − b − a 00 = 0 ) + 1 ) 0 + b + ( a + ( b − 1 − a − b 00 = 0 ) + 1 ) 2 12 ∆ y 2 12 ∆ x 1 ∆ x 1 ∆ y a + ( b + 0 − a − 2 − b − = 0 + 2 ) b + ( a + a 10 0 − b − 2 − a − b 01 = 0 + 2 ) 30 ∆ y 30 ∆ x 1 ∆ x 3 ∆ x 1 ∆ y 3 ∆ y ( b + ( b + ( a + ( a + 1 − b − 3 − b − 1 − a − 3 − a − a 20 = − 1 ) + 3 ) b 02 = − 1 ) + 3 ) 2 ∆ y 140 ∆ y 2 ∆ x 140 ∆ x 1 1 ∆ x 3 + b + ( b + ( b − 2 − b − b 30 = 3 ) a 30 = − 2 ) 2 3 ∆ y 1 ∆ y 1 ( a + 2 − a − 3 + a + b 03 = − 2 ) ( a − a 03 = 3 ) 3 ∆ x 2 b + a + 2 − b − 2 − a − b 21 = = a 12 2 2 b + 3 − b − 1 ∆ x b 31 = ( b + 3 = − 3 − b − 3 ) a 40 4 ∆ y 1 ∆ y ( a + 3 − a − b 04 = − 3 ) a + 3 − a − 4 ∆ x a 13 = 3 18 / 38
Recommend
More recommend
