Exactly divergence-free methods for induction and related equations - PowerPoint PPT Presentation

Exactly divergence-free methods 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 7 May 2018 Supported by Airbus Foundation Chair at TIFR-CAM, Bangalore http://math.tifrbng.res.in/airbus-chair 1 / 41

Outline 1 DG scheme using Raviart-Thomas polynomials 2 Divergence-free reconstruction approach Joint work with Dinshaw S. Balsara, Univ. of Notre Dame Rakesh Kumar, TIFR-CAM Arijit Hazra, TIFR-CAM 2 / 41

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

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

I. DG scheme based on Raviart-Thomas polynomials 5 / 41

Approximation of magnetic field Find B ∈ H (div , Ω) H (div , Ω) = { B ∈ L 2 (Ω) : div( B ) ∈ L 2 (Ω) } To approximate functions in H (div , Ω) on a mesh T h , the normal trace B · n must be continuous. P k ( x ) , P k ( y ) : 1-D polynomials of degree at most k Q r,s ( x, y ) : tensor product polynomials 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) 6 / 41

Approximation of magnetic field Two consequences: • For any B h ∈ RT k , we have div( B h ) ∈ Q k,k ( x, y ) =: Q k ( x, y ) • The restriction of B h · n 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 7 / 41

Approximation of magnetic field 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 B x ∈ Q 1 , 0 B y ∈ Q 0 , 1 Location of dofs of Raviart-Thomas polynomial for k = 0 8 / 41

Approximation of magnetic field B x ∈ Q 2 , 1 B y ∈ Q 1 , 2 Location of dofs of Raviart-Thomas polynomial for k = 1 B x ∈ Q 3 , 2 B y ∈ Q 2 , 3 Location of dofs of Raviart-Thomas polynomial for k = 2 9 / 41

Approximation of magnetic field 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 [ ? ], [ ? ] � 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 (Ω) 10 / 41

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 11 / 41

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. 12 / 41

Theorem If all the moments are zero for any cell C , then B h ≡ 0 inside that cell. Theorem Let B h ∈ RT k satisfy the moments � � B h x φ d y = B x φ d y ∀ φ ∈ P k ( y ) (1) e ∓ e ∓ x x � � B h y φ d x = B y φ d x ∀ φ ∈ P k ( x ) (2) e ∓ e ∓ y y � � B h x ψ d x d y = B x ψ d x d y ∀ ψ ∈ ∂ x Q k ( x, y ) (3) C C � � B h y ψ d x d y = B y ψ d x d y ∀ ψ ∈ ∂ y Q k ( x, y ) (4) C C for a given vector field B ∈ H ( div , Ω) . If div ( B ) ≡ 0 then div ( B h ) ≡ 0 . 13 / 41

Model problem ∂ B ∂t + ∇ × E = − M Divergence evolves according to ∂ ∂t (div B ) + ∇ · M = 0 (5) If M ≡ 0 , then div( B ) = constant. In 2-D ∂B x ∂t + ∂E z ∂B y ∂t − ∂E z ∂y = − M x , ∂x = − M y and for the induction equation E z = v y B x − v x B y 14 / 41

DG scheme for the induction equation Constructing B h from the edge and cell moments allowed us to get divergence-free approximation. We will construct a scheme to evolve the same moments in time . e + Consider the right edge e + 2 y 3 x . Multiply by test function φ ∈ P k ( y ) � ∂B x � � ∂t + ∂E � e + e − φ d y = − M x φ d y C x x ∂y e + e + x x and integrate by parts in second term 0 1 e − y � ∂B x � E ∂φ � ∂t φ d y − ∂y d y + ( Eφ ) 3 − ( Eφ ) 1 = − M x φ d y e + e + e + x x x 15 / 41

DG scheme for the induction equation Edge moments are evolved by ∂B h E ∂φ � � � ˆ ∂y d y + [ ˜ ˆ x ∂t φ d y − Eφ ] e ∓ x = − M x φ d y, ∀ φ ∈ P k ( y ) e ∓ e ∓ e ∓ x x x ∂B h � � E ∂φ � y ˆ ∂x d x − [ ˜ ˆ ∂t φ d x + Eφ ] e ∓ y = − M y φ d y, ∀ φ ∈ P k ( x ) e ∓ e ∓ e ∓ y y y where ˆ E = numerical flux from a 1-D Riemann solver required on the faces ˜ E = numerical flux from a multi-D Riemann solver needed at vertices [ ˜ x = ( ˜ Eφ ) 2 − ( ˜ [ ˜ x = ( ˜ Eφ ) 3 − ( ˜ Eφ ] e − Eφ ) 0 , Eφ ] e + Eφ ) 1 [ ˜ y = ( ˜ Eφ ) 1 − ( ˜ [ ˜ y = ( ˜ Eφ ) 3 − ( ˜ Eφ ] e − Eφ ) 0 , Eφ ] e + Eφ ) 2 16 / 41

DG scheme for the induction equation The cells moments are evolved by the following standard DG scheme ∂B h � � E ∂ψ � � x ˆ ∂t ψ d x d y − ∂y d x d y + Eψn y d s = − M x ψ d x d y, ∀ ψ ∈ ∂ x Q k ( x, y ) C C ∂C C ∂B h � � E ∂ψ � � y ˆ ∂t ψ d x d y + ∂x d x d y − Eψn x d s = − M y ψ d x d y, ∀ ψ ∈ ∂ y Q k ( x, y ) C C ∂C C Note that the same 1-D numerical flux ˆ E is used in both the edge and cell moment equations whereas the vertex numerical flux ˜ E is needed only in the edge moment equations. 17 / 41

Theorem If M = 0 , the DG scheme preserves the divergence of the magnetic field. Proof: For any φ ∈ Q k take ψ = ∂ x φ and ψ = ∂ y φ in the two cell moment equations respectively and add them together to obtain � ∂t ∂ x φ + ∂B h � ∂B h � � y x ˆ ∂t ∂ y φ d x d y − E ( n x ∂ y φ − n y ∂ x φ ) d s = 0 C ∂C Two of the cell integrals cancel since ∂ x ∂ y φ = ∂ y ∂ x φ . Performing an integration by parts in the first term, we obtain � φ ∂ � φ ∂ � ˆ ∂t div( B h ) d x d y + ∂t ( B h · n ) d s − − E ( n x ∂ y φ − n y ∂ x φ ) d s = 0 C ∂C ∂C The last two terms can be re-arranged as φ∂B h φ∂B h φ∂B h φ∂B h � � � � x x y y ∂t d y − ∂t d y + ∂t d x − ∂t d x e + e + e − e − x x y y � � � � ˆ ˆ ˆ ˆ − E∂ y φ d y + E∂ y φ d y + E∂ x φ d x − E∂ x φ d x e + e + e − e − x x y y 18 / 41

The restriction of φ on each edge is a one dimensional polynomial of degree k . We can use the edge moment equations to simplify as = − [ ˜ x + [ ˜ x + [ ˜ y − [ ˜ Eφ ] e + Eφ ] e − Eφ ] e + Eφ ] e − y = − ( ˜ Eφ ) 3 + ( ˜ Eφ ) 1 + ( ˜ Eφ ) 2 − ( ˜ Eφ ) 0 + ( ˜ Eφ ) 3 − ( ˜ Eφ ) 2 − ( ˜ Eφ ) 1 + ( ˜ Eφ ) 0 = 0 Hence we have � φ ∂ ∂t div( B h ) d x d y = 0 ∀ φ ∈ Q k ( x, y ) C Since div( B h ) ∈ Q k ( x, y ) we conclude that the divergence is preserved by the numerical scheme. 19 / 41

Theorem ( M � = 0 ) The divergence evolves consistently with equation (5) in the sense that � φ ∂ � � ∂t div ( B h ) d x d y − φ ˆ M ·∇ φ d x d y + M · n d s = 0 , ∀ φ ∈ Q k C C ∂C Since div ( B h ) ∈ Q k , we can expect div ( B h ) to be accurate to O ( h k +1 ) . In case of Maxwell equations, if we want to compute charge density ρ = ∇ · D we can get ρ to same accuracy as the vector fields. 20 / 41

Numerical fluxes Using the zero divergence condition, rewrite the induction equation ∂B x ∂v y ∂v x ∂B y ∂v x ∂v y ∂t + v ·∇ B x + B x ∂y − B y ∂y = 0 , ∂t + v ·∇ B y + B y ∂x − B x ∂x = 0 The characteristics are the integral curves of v . The 1-D numerical flux is given by � E L if v · n > 0 ˆ E = E R otherwise For example, across the face e ∓ x , the flux is given by � v y B x − v x B L if v x > 0 ˆ y E = v y B x − v x B R otherwise y 21 / 41