Well-balanced DG scheme for Euler equations with gravity Praveen - PowerPoint PPT Presentation

Well-balanced DG scheme for Euler equations with gravity Praveen Chandrashekar praveen@tifrbng.res.in Center for Applicable Mathematics Tata Institute of Fundamental Research Bangalore 560065 Higher Order Numerical Methods for Evolutionary PDEs: Applied Mathematics Meets Astrophysical Applications (15w5134) BIRS, Banff Centre, 10–15 May, 2015 Supported by Airbus Foundation Chair at TIFR-CAM, Bangalore 1 / 55

Euler equations with gravity   0 ∂ q ∂t + ∂ f  ∂ Φ ∂x = − ρ  ∂x ρu where     ρ ρu  ,  , p + ρu 2 q = ρu f = Φ = gravitational potential   E ( E + p ) u Calorically ideal gas � � E − 1 γ = c p 2 ρu 2 p = ( γ − 1) , > 1 c v 2 / 55

Hydrostatic solutions • Fluid at rest u e = 0 • Mass and energy equation satisfied ∂ t ρ + ∂ x ( ρu ) = 0 , ∂ t E + ∂ x [( E + p ) u ] = − ρu∂ x Φ • Momentum equation d p e dΦ d x = − ρ e (1) d x • Need additional assumptions to solve this equation 3 / 55

Hydrostatic solutions: Isothermal case • Thermally ideal gas model p = ρRT, R = gas constant • Isothermal hydrostatic state, i.e., T e ( x ) = T e = const , then � Φ( x ) � p e ( x ) exp = const (2) RT e Density ρ e ( x ) = p e ( x ) RT e 4 / 55

Hydrostatic solutions: Polytropic case • Polytropic hydrostatic state p e ρ − ν = α = const , ν > 1 constant (3) e From (1) and (3), we obtain ναρ ν − 1 ( x ) e + Φ( x ) = β = const (4) ν − 1 5 / 55

Scope of present work • Nodal DG scheme ◮ Gauss-Lobatto-Legendre nodes ◮ Arbitrary quadrilateral cells in 2-D • Well-balanced for hydrostatic solutions ◮ isothermal solutions under ideal gas ◮ polytropic solutions (with Markus Zenk) • Any consistent numerical flux ◮ Discontinuous density: contact preserving flux 6 / 55

Finite element method Conservation law with source term stationary solution q e q t + f ( q ) x = s ( q ) f ( q e ) x = s ( q e ) Weak formulation : Find q ( t ) ∈ V such that d d t ( q, ϕ ) + a ( q, ϕ ) = ( s ( q ) , ϕ ) , ∀ ϕ ∈ V FEM with quadrature : Find q h ( t ) ∈ V h such that d d t ( q h , ϕ h ) h + a h ( q h , ϕ h ) = ( s h ( q h ) , ϕ h ) h , ∀ ϕ h ∈ V h 7 / 55

Finite element and well-balanced Stationary solution is not a polynomial, q e / ∈ V h . Let q e,h = Π h ( q e ) , Π h : V → V h , interpolation or projection FEM is well-balanced if a h ( q e,h , ϕ h ) = ( s h ( q e,h ) , ϕ h ) h , ∀ ϕ h ∈ V h because if q h (0) = q e,h = ⇒ q h ( t ) = q e,h ∀ t 8 / 55

Mesh and basis functions • Partition domain into disjoint cells C i = ( x i − 1 2 , x i + 1 2 ) , ∆ x i = x i + 1 2 − x i − 1 2 • Approximate solution inside each cell by a polynomial of degree N I i − 1 I i I i − 1 9 / 55

Mesh and basis functions • Map C i to a reference cell, say [0 , 1] x = ξ ∆ x i + x i − 1 (5) 2 • On reference cell, ξ j , 0 ≤ j ≤ N are Gauss-Lobatto-Legendre nodes, roots of (1 − ξ 2 ) P ′ N ( ξ ) in [ − 1 , +1] • ℓ j ( ξ ) = Lagrange polynomials using GLL points ℓ j ( ξ k ) = δ jk , 0 ≤ j, k ≤ N • Basis functions in physical coordinates ϕ j ( x ) = ℓ j ( ξ ) , 0 ≤ j ≤ N • Derivative of ϕ j : apply the chain rule of differentiation d j ( ξ ) d ξ 1 d xϕ j ( x ) = ℓ ′ ℓ ′ d x = j ( ξ ) ∆ x i 10 / 55

Mesh and basis functions • x j ∈ C i denote the physical locations of the GLL points x j = ξ j ∆ x i + x i − 1 2 , 0 ≤ j ≤ N 11 / 55

Discontinuous Galerkin Scheme Consider the single conservation law with source term ∂q ∂t + ∂f ∂x = s ( q, x ) Solution inside cell C i is polynomial of degree N N � q h ( x, t ) = q j ( t ) ϕ j ( x ) , q j ( t ) = q h ( x j , t ) j =0 Approximate the flux N N � � f ( q h ) ≈ f h ( x, t ) = f ( q h ( x j , t )) ϕ j ( x ) = f j ( t ) ϕ j ( x ) j =0 j =0 12 / 55

Discontinuous Galerkin Scheme Gauss-Lobatto-Legendre quadrature N � � φ ( x ) ψ ( x )d x ≈ ( φ, ψ ) h = ∆ x i ω q φ ( x q ) ψ ( x q ) C i q =0 Semi-discrete DG: For 0 ≤ j ≤ N d d t ( q h , ϕ j ) h + ( ∂ x f h , ϕ j ) h +[ ˆ 2 − f h ( x − 2 )] ϕ j ( x − f i + 1 2 ) i + 1 i + 1 (6) − [ ˆ 2 − f h ( x + 2 )] ϕ j ( x + f i − 1 2 ) = ( s h , ϕ j ) h i − 1 i − 1 where ˆ 2 = ˆ f ( q − 2 , q + f i + 1 2 ) is a numerical flux function. This is also i + 1 i + 1 called a DG Spectral Element Method . 13 / 55

Numerical flux Consistency of numerical flux ˆ f ( q , q ) = f ( q ) Def: Contact property The numerical flux ˆ f is said to satisfy contact property if for any two states w L = [ ρ L , 0 , p ] w R = [ ρ R , 0 , p ] and we have ˆ f ( q ( w L ) , q ( w R )) = [0 , p, 0] ⊤ • states w L , w R in the above definition correspond to a stationary contact discontinuity. • Contact Property = ⇒ exactly supports stationary contact • Examples: Roe, HLLC, etc. 14 / 55

Approximation of source term: isothermal case Let ¯ T i = temperature corresponding to the cell average value in cell C i Rewrite the source term in the momentum equation as (Xing & Shu) � Φ � ∂ � � s = − ρ∂ Φ − Φ ∂x = ρR ¯ T i exp ∂x exp R ¯ R ¯ T i T i Source term approximation: For x ∈ C i � ∂ � Φ( x ) N � � − Φ( x j ) � s h ( x ) = ρ h ( x ) R ¯ T i exp exp ϕ j ( x ) (7) R ¯ R ¯ T i ∂x T i j =0 Source term in the energy equation 1 ( ρu ) h s h ρ h 15 / 55

Approximation of source term: polytropic case Define H ( x ) inside each cell C i as � ν − 1 � ν H ( x ) = ν − 1 ln ( β i − Φ( x )) , x ∈ C i να i α i , β i : constants to be chosen. Rewrite source term s ( x ) = − ρ∂ Φ ∂x = ν − 1 ρ ( β i − Φ( x )) exp( − H ( x )) ∂ ∂x exp( H ( x )) , x ∈ C i ν The source term is approximated as N s h ( x ) = ν − 1 ρ h ( x )( β i − Φ( x )) exp( − H ( x )) ∂ � exp( H ( x j )) ϕ j ( x ) ν ∂x j =0 � � ν p j α i = p j ∗ ρ − ν β i = max + Φ( x j ) , j ∗ ν − 1 ρ j 0 ≤ j ≤ N 16 / 55

Well-balanced property Well-balanced property Let the initial condition to the DG scheme (6), (7) be obtained by interpolating the isothermal/polytropic hydrostatic solution corresponding to a continuous gravitational potential Φ . Then the scheme (6), (7) preserves the initial condition under any time integration scheme. Proof : For continuous hydrostatic solution, q h (0) is continuous. By flux consistency ˆ ˆ 2 − f h ( x − 2 − f h ( x + f i + 1 2 ) = 0 , f i − 1 2 ) = 0 i + 1 i − 1 Above is true even if density is discontinuous, provided flux satisfies contact property. = ⇒ density and energy equations are well-balanced 17 / 55

Well-balanced property Momentum eqn: flux f h has the form N � f h ( x, t ) = p j ( t ) ϕ j ( x ) , p j = pressure at the GLL point x j j =0 Isothermal initial condition, ¯ T i = T e = const . The source term evaluated 18 / 55

Well-balanced property at any GLL node x k is given by � ∂ � N � Φ( x k ) � − Φ( x j ) � s h ( x k ) = ρ h ( x k ) RT e exp exp ∂xϕ j ( x k ) RT e RT e � �� � j =0 p k � ∂ � N � Φ( x k ) � − Φ( x j ) � = p k exp exp ∂xϕ j ( x k ) RT e RT e j =0 � ∂ N � Φ( x k ) � � − Φ( x j ) � = p k exp exp ∂xϕ j ( x k ) RT e RT e j =0 � ∂ N � Φ( x j ) � � − Φ( x j ) � = p j exp exp ∂xϕ j ( x k ) RT e RT e j =0 N ∂xϕ j ( x k ) = ∂ ∂ � = p j ∂xf h ( x k ) j =0 19 / 55

Well-balanced property Since ∂ x f h ( x k ) = s h ( x k ) at all the GLL nodes x k we can conclude that ( ∂ x f h , ϕ j ) h = ( s h , ϕ j ) h , 0 ≤ j ≤ N = ⇒ scheme is well-balanced for the momentum equation The proof for polytropic case is similar. 20 / 55

2-D Euler equations with gravity ∂ q ∂t + ∂ f ∂x + ∂ g ∂y = s q = vector of conserved variables, ( f , g ) = flux vector and s = source term, given by       ρ ρu ρv p + ρu 2 ρu ρuv       q =  , f =  , g =       p + ρv 2 ρv ρuv     E ( E + p ) u ( E + p ) v   0 − ρ ∂ Φ   ∂x   s = − ρ ∂ Φ   ∂y   � � u ∂ Φ ∂x + v ∂ Φ − ∂y 21 / 55

2-D hydrostatic solution: Isothermal case Momentum equation ∇ p e = − ρ e ∇ Φ Assuming the ideal gas equation of state p = ρRT and a constant temperature T = T e = const , we get � Φ( x, y ) � p e ( x, y ) exp = const (8) RT e We will exploit the above property of the hydrostatic state to construct the well-balanced scheme. 22 / 55

Mesh and basis functions A quadrilateral cell K and reference cell ˆ K = [0 , 1] × [0 , 1] y η 4 3 K 4 3 T K 2 1 ˆ K x ξ 1 2 1-D GLL points ξ r ∈ [0 , 1] , 0 ≤ r ≤ N 23 / 55

Mesh and basis functions Tensor product of GLL points N = 1 N = 2 N = 3 Basis functions in cell K : i ↔ ( r, s ) ϕ K i ( x, y ) = ℓ r ( ξ ) ℓ s ( η ) , ( x, y ) = T K ( ξ, η ) Computation of ∂ x f h requires derivatives of the basis functions ∂ r ( ξ ) ℓ s ( η ) ∂ξ s ( η ) ∂η ∂xϕ K i ( x, y ) = ℓ ′ ∂x + ℓ r ( ξ ) ℓ ′ ∂x 24 / 55