Divergence-free discontinuous Galerkin method for ideal compressible - PowerPoint PPT Presentation

Divergence-free discontinuous Galerkin method for ideal compressible MHD equations Praveen Chandrashekar praveen@math.tifrbng.res.in http://cpraveen.github.io Center for Applicable Mathematics Tata Institute of Fundamental Research Bangalore-560065, India http://math.tifrbng.res.in Oberseminar, Dept. of Mathematics, University of W¨ urzburg 4 June 2019 1 / 35

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 / 35

Ideal compressible MHD equations Nonlinear hyperbolic system Compressible Euler equations with Lorentz force ∂ρ ∂t + ∇ · ( ρ v ) = 0 ∂ ( ρ v ) + ∇ · ( P I + ρ v ⊗ v − B ⊗ B ) = 0 ∂t ∂E ∂t + ∇ · (( E + P ) v + ( v · B ) B ) = 0 ∂ B ∂t − ∇ × ( v × B ) = 0 Magnetic monopoles do not exist: = ⇒ ∇ · B = 0 3 / 35

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 Wu) 4 / 35

Objectives • Based on conservation form of the equations • Upwind-type schemes using Riemann solvers • Divergence-free schemes for Maxwell’s and compressible MHD ◮ Cartesian grids at present ◮ Divergence preserving schemes (RT) ◮ Divergence-free reconstruction (BDM) • High order accurate ◮ discontinuous-Galerkin • Non-oscillatory schemes for MHD ◮ using limiters • Explicit time stepping • Based on previous work for induction equation ◮ J. Sci. Comp., Vol. 79, pp, 79-102, 2019 5 / 35

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 / 35

MHD equations in 2-D ∂ U ∂t + ∂ F x ∂x + ∂ F y ∂y = 0       ρ ρv x ρv y P + ρv 2 x − B 2 ρv x v y − B x B y ρv x x            P + ρv 2 y − B 2  ρv x v y − B x B y ρv y y             ρv x v z − B x B z ρv y v z − B y B z ρv z       U = , F x = , F y =       E ( E + P ) v x − B x ( v · B ) ( E + P ) v y − B y ( v · B )             B x 0 E z             B y − E z 0       B z v x B z − v z B x v y B z − v z B y where P = p + 1 γ − 1 + 1 p 2 ρ | v | 2 + 1 2 | B | 2 , 2 | B | 2 B = ( B x , B y , B z ) , E = E z is the electric field in the z direction E z = − ( v × B ) z = v y B x − v x B y 7 / 35

MHD equations in 2-D Split into two parts U = [ ρ, ρ v , E , B z ] ⊤ , B = ( B x , B y ) ∂ U ∂B x ∂t + ∂E z ∂B y ∂t − ∂E z ∂t + ∇ · F ( U , B ) = 0 , ∂y = 0 , ∂x = 0 The fluxes F = ( F x , F y ) are of the form     ρv x ρv y P + ρv 2 x − B 2 ρv x v y − B x B y     x     P + ρv 2 y − B 2 ρv x v y − B x B y     y F x = , F y =     ρv x v z − B x B z ρv y v z − B y B z         ( E + P ) v x − B x ( v · B ) ( E + P ) v y − B y ( v · B )     v x B z − v z B x v y B z − v z B y 8 / 35

Approximation of magnetic field If we want ∇ · 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 9 / 35

Approximation spaces: Degree k ≥ 0 Map cell K to reference cell ˆ K = [ − 1 2 , + 1 2 ] × [ − 1 2 , + 1 2 ] P r ( ξ ) = span { 1 , ξ, ξ 2 , . . . , ξ r } , Q r,s ( ξ, η ) = P r ( ξ ) ⊗ P s ( η ) Hydrodynamic variables in each cell k k � � U ( ξ, η ) = U ij φ i ( ξ ) φ j ( η ) ∈ Q k,k i =0 j =0 Normal component of B on faces k � on vertical faces : b x ( η ) = a j φ j ( η ) ∈ P k ( η ) j =0 k � on horizontal faces : b y ( ξ ) = b j φ j ( ξ ) ∈ P k ( ξ ) j =0 { φ i ( ξ ) } are orthogonal polynomials on [ − 1 2 , + 1 2 ] , with degree φ i = i . 10 / 35

Approximation spaces: Degree k ≥ 0 For k ≥ 1 ,define certain cell moments � + 1 � + 1 1 2 2 α ij = α ij ( B x ) := B x ( ξ, η ) φ i ( ξ ) φ j ( η ) d ξ d η, 0 ≤ i ≤ k − 1 , 0 ≤ j ≤ k m ij − 1 − 1 2 2 � + 1 � + 1 1 2 2 β ij = β ij ( B y ) := B y ( ξ, η ) φ i ( ξ ) φ j ( η ) d ξ d η, 0 ≤ i ≤ k, 0 ≤ j ≤ k − 1 m ij − 1 − 1 2 2 � + 1 � + 1 � + 1 2 2 2 [ φ i ( ξ ) φ j ( η )] 2 d ξ d η = m i m j , [ φ i ( ξ )] 2 d ξ m ij = m i = − 1 − 1 − 1 2 2 2 α 00 , β 00 are cell averages of B x , B y Solution variables { U ( ξ, η ) } , { b x ( η ) } , { b y ( ξ ) } , { α, β } The set b x , b y , α, β are the dofs for the Raviart-Thomas space. 11 / 35

RT reconstruction: b ± x ( η ) , b ± y ( ξ ) , α, β → B ( ξ, η ) Given b ± x ( η ) ∈ P k and b ± 2 3 y ( ξ ) ∈ P k , b + y ( ξ ) and set of cell moments x ( η ) α x ( η ) { α ij , 0 ≤ i ≤ k − 1 , 0 ≤ j ≤ k } b − b + β { β ij , 0 ≤ i ≤ k, 0 ≤ j ≤ k − 1 } b − y ( ξ ) 0 1 Find B x ∈ Q k +1 ,k and B y ∈ Q k,k +1 such that 2 , η ) = b ± 2 ) = b ± B x ( ± 1 η ∈ [ − 1 2 , 1 B y ( ξ, ± 1 ξ ∈ [ − 1 2 , 1 x ( η ) , 2 ] , y ( ξ ) , 2 ] � + 1 � + 1 1 2 2 B x ( ξ, η ) φ i ( ξ ) φ j ( η ) d ξ d η = α ij , 0 ≤ i ≤ k − 1 , 0 ≤ j ≤ k m ij − 1 − 1 2 2 � + 1 � + 1 1 2 2 B y ( ξ, η ) φ i ( ξ ) φ j ( η ) d ξ d η = β ij , 0 ≤ i ≤ k, 0 ≤ j ≤ k − 1 m ij − 1 − 1 2 2 (1) ∃ unique solution. (2) Data div-free = ⇒ reconstructed B is div-free. 12 / 35

DG scheme for B on faces On every vertical face of the mesh � + 1 � + 1 ∂t φ i d η − 1 d φ i d η d η + 1 ∂b x 2 2 ˆ ∆ y [ ˜ E z E z φ i ] = 0 , 0 ≤ i ≤ k ∆ y − 1 − 1 2 2 On every horizontal face of the mesh � + 1 � + 1 ∂b y ∂t φ i d ξ + 1 d φ i d ξ d ξ − 1 2 2 ˆ ∆ x [ ˜ E z E z φ i ] = 0 , 0 ≤ i ≤ k ∆ x − 1 − 1 2 2 ˆ E z : on face, 1-D Riemann solver ˜ E z : at vertex, 2-D Riemann solver 13 / 35

DG scheme for B on cells � + 1 � + 1 m ij d α ij ∂B x 2 2 = ∂t φ i ( ξ ) φ j ( η ) d ξ d η d t − 1 − 1 2 2 � + 1 � + 1 − 1 ∂E z 2 2 = ∂η φ i ( ξ ) φ j ( η ) d ξ d η ∆ y − 1 − 1 2 2 � + 1 − 1 2 [ ˆ 2 ) − ˆ E z ( ξ, 1 2 ) φ i ( ξ ) φ j ( 1 E z ( ξ, − 1 2 ) φ i ( ξ ) φ j ( − 1 = 2 )] d ξ ∆ y − 1 2 � + 1 � + 1 + 1 2 2 E z ( ξ, η ) φ i ( ξ ) φ ′ j ( η ) d ξ d η, 0 ≤ i ≤ k − 1 , 0 ≤ j ≤ k ∆ y − 1 − 1 2 2 Not a Galerkin method, test functions ( Q k − 1 ,k ) different from trial functions ( Q k +1 ,k ) 14 / 35

DG scheme for U on cells For each test function Φ( ξ, η ) = φ i ( ξ ) φ j ( η ) ∈ Q k,k � 1 � + 1 � + 1 � + 1 � + 1 ∂ U c � ∆ x F x ∂ Φ ∆ y F y ∂ Φ 1 2 2 2 2 ∂t Φ( ξ, η ) d ξ d η − ∂ξ + d ξ d η ∂η − 1 − 1 − 1 − 1 2 2 2 2 � + 1 � + 1 + 1 1 2 2 F + ˆ F − ˆ x Φ( 1 x Φ( − 1 2 , η ) d η − 2 , η ) d η ∆ x ∆ x − 1 − 1 2 2 � + 1 � + 1 + 1 1 2 2 F + ˆ ˆ y Φ( ξ, 1 F − y Φ( ξ, − 1 2 ) d ξ − 2 ) d ξ = 0 ∆ y ∆ y − 1 − 1 2 2 15 / 35

DG scheme for U on cells U n B n x B n y b + y U w U c U e B w B c B e b − b + x x x x x B w B c B e y y y b − y U s B s x B s y F x = F x ( U c , B c x , B c F y = F y ( U c , B c x , B c y ) , y ) F + ˆ x = ˆ F x (( U c , b + y ) , ( U e , b + ˆ x = ˆ x , B c x , B e F − F x (( U w , b − x , B w y ) , ( U c , b − x , B c y )) , y )) F + ˆ y = ˆ F y (( U c , B c x , b + y ) , ( U n , B n x , b + F − ˆ y = ˆ F y (( U s , B s x , b − y ) , ( U c , B c x , b − y )) , y )) 16 / 35

Constraints on B Definition (Strongly divergence-free) We will say that a vector field B defined on a mesh is strongly divergence-free if 1 ∇ · B = 0 in each cell K ∈ T h 2 B · n is continuous at each face F ∈ T h Theorem (1) The DG scheme satisfies � d ( ∇ · B ) φ d x d y = 0 , ∀ φ ∈ Q k,k d t K and since ∇ · B ∈ Q k,k = ⇒ ∇ · B = constant wrt time. (2) If ∇ · B ≡ 0 at t = 0 = ⇒ ∇ · B ≡ 0 for t > 0 17 / 35