An Operator Splitting Based Stochastic Galerkin Method for Nonlinear - PowerPoint PPT Presentation

An Operator Splitting Based Stochastic Galerkin Method for Nonlinear Systems of Hyperbolic Conservation Laws with Uncertainty Alexander Kurganov Tulane University, Mathematics Department www.math.tulane.edu/ ∼ kurganov Supported by NSF

joint work with Alina Chertock, North Carolina State University, USA Shi Jin, University of Wisconsin – Madison, USA 2

Conservation/Balance Laws with Uncertainties x ∈ R , t > 0 , z ∈ Ω ⊂ R d U t + F ( U , x, z ) x = R ( U , x, z ) , U = U ( x, t, z ) is the unknown vector function x : spatial variable t : time variable z : random variable F : flux vector function R : source term Uncertainties can appear in the source terms, equations of state, initial or boundary data due to empirical approximations or measuring errors

Quantifying Uncertainties – gPC Approach Polynomial chaos or generalized polynomial chaos (gPC) approach: • Non-intrusive gPC method – solves the original problem at selected sampling points, thus one can use the deterministic code, and then use interpolation and quadrature rules to numerically evaluate the statistical moments [Xiu, Hesthaven; 2005] [Mishra, Schwab, Sukys; 2012] • Intrusive gPC method – uses the Galerkin approximation, which results in a system of deterministic equations, solving which will give the stochastic moments of the solution of the original uncertain problem 4

– Pros: lower computational cost; theoretical advantages; [Elman, Miller, Phipps, Tuminaro; 2011] – Cons: extra efforts are needed in order to obtain well-behaved discrete systems [Xui; 2010] [Tryoen, Le Maitre, Ndjinga, Ern; 2010] [Despr´ es, Po¨ ette, Lucor; 2013] [Pettersson, Iaccarino, Nordstr¨ om; 2014, 2015] [Hu, Jin, Xiu; 2015] 5

The gPC-SG Method – An Overview x ∈ R , t > 0 , z ∈ Ω ⊂ R d U t + F ( U , x, z ) x = R ( U , x, z ) , The solution is sought in terms of an orthogonal polynomial series in z M − 1 � d + N � ˆ � U ( x , t, z ) ≈ U N ( x , t, z ) = U i ( x , t )Φ i ( z ) , M = d i =0 • { Φ i ( z ) } are multidimensional polynomials of degree up to N of z : � � P d � Φ i ( z )Φ ℓ ( z ) µ ( z ) d z = δ iℓ , 0 ≤ i, ℓ ≤ M − 1 M = dim N Ω • µ ( z ): probability density function of z • δ iℓ : Kronecker symbol • The choice of the orthogonal polynomials depends on the distribution function of z . For example: – a Gaussian distribution defines the Hermite polynomials – a uniform distribution defines the Legendre polynomials 6

The gPC-SG method seeks to satisfy the system in a weak form by ensuring that the residual is orthogonal to the gPC polynomial space. M − 1 ˆ Substituting U N ( x , t, z ) = U i ( x , t )Φ i ( z ) � i =0 into the governing system U t + F ( U , x, z ) x = R ( U , x, z ) and using the Galerkin projection yield ( ˆ U i ) t + ( ˆ F i ) x = ˆ R i , 0 ≤ i ≤ M − 1 where � M − 1 � � ˆ ˆ � F i = U j ( x , t )Φ j ( z ) , x, z Φ i ( z ) µ ( z ) d z F j =0 Ω � M − 1 � � ˆ ˆ � R i = U i ( x , t )Φ i ( z ) , x, z Φ i ( z ) µ ( z ) d z R j =0 Ω 7

The gPC-SG Method – Challenges ( ˆ U i ) t + ( ˆ F i ) x = ˆ U t + F ( U , x, z ) x = R ( U , x, z ) 0 ≤ i ≤ M − 1 R i Linear Hyperbolic Hyperbolic Nonlinear Symmetric Hyperbolic Nonlinear Nonsymmetric ? • Our goal: Introduce an operator splitting for the original hyperbolic system, which will guarantee that the gPC-SG discretization of each of the split subsystems always results in a globally hyperbolic system • Our strategy: generic, but the splitting is problem specific • Our examples: the compressible Euler equations and the shallow water equations 8

1-D Compressible Euler Equations  ρ t + m x = 0   m t + ( ρu 2 + p ) x = 0    E t + ( u ( E + p )) x = 0  • ρ : density • u : velocity, m = ρu : momentum • E : total energy E − 1 � 2 ρu 2 � • p : pressure with the equation of state p = ( γ − 1) • γ : specific heat ratio We assume here that the data may depend on random variable z , i.e., ρ ( x, 0 , z ) = ρ 0 ( x, z ) , u ( x, 0 , z ) = u 0 ( x, z ) , p ( x, 0 , z ) = p 0 ( x, z ) , γ = γ ( z ) Uncertainty may also arise from boundary data and other terms

1-D Euler Equations – Numerical Challenges  ρ t + m x = 0   m t + ( ρu 2 + p ) x = 0 �  λ = u, u ± c, c = γp/ρ   E t + ( u ( E + p )) x = 0  A direct application of the gPC-SG method to the system may fail due to the loss of hyperbolicity after the gPC-SG discretization Operator Splitting : • Linear hyperbolic system • Two degenerate nonlinear hyperbolic systems which are effectively scalar equations The gPC-SG approximation is guaranteed to maintain the hyperbolicity for each of the subsystems 10

1-D Euler Equations – Operator Splitting  ρ t + m x = 0    ( I ) m t + (( γ − 1) E + am ) x = 0   E t − ( aE ) x = 0   ρ t = 0    · m 2  � � 3 − γ   ( II ) m t + ρ − am = 0 2 x     E t = 0    ρ t = 0     m t = 0   ( III ) · m 2 � � � � m γE − γ − 1   E t + + aE = 0   ρ 2 ρ   x • We choose: −| a | ≤ u − c < u + c ≤ | a | : subcharacteristic condition a = ± sup(max {| u | + c, γu, (3 − γ ) u } ) : convection coefficient should not change sign 11

Strang Splitting U t + F I ( U ) x = 0 → S I U t + F II ( U ) x = 0 → S II U t + F III ( U ) x = 0 → S III Here   0     r m · m 2 3 − γ   U = F I = ( γ − 1) E + am F II = m  ,  , ρ − am       2     E − aE 0   0 0   F III =   � �  · m 2  γE − γ − 1 m + aE   ρ 2 ρ • Assume that the solution of the original system is available at time t • Introduce a (small) time step ∆ t • One time step of the second-order Strang splitting method: U ( x , t + ∆ t, z ) = S I (∆ t/ 2) S II (∆ t/ 2) S III (∆ t ) S II (∆ t/ 2) S I (∆ t/ 2) U ( x , t, z )

Operator Splitting – Numerical Validation • We consider the Sod shock tube problem – pure deterministic problem: � 1 , � 1 , x < 0 . 5 , x < 0 . 5 ρ 0 ( x ) = u 0 ( x ) ≡ 0 , p 0 ( x ) = 0 . 125 , x > 0 . 5 , 0 . 1 , x > 0 . 5 • We run numerical simulations for both the unsplit and split systems • We compare the results computed by the central-upwind scheme – computational domain [0,1] – non-reflecting boundary conditions – uniform grid with ∆ x = 1 / 400 – final time t = 0 . 1644 13

ρ (top left), m (top right) and E (bottom) 14

1-D Euler Equations – The gPC-SG Approximation  ρ t = 0   ρ t + m x = 0   · m 2   � � 3 − γ     m t + (( γ − 1) E + am ) x = 0 m t + ρ − am = 0 2 x     E t − ( aE ) x = 0    E t = 0 ,    ρ t = 0     m t = 0   · m 2 � � � � γE − γ − 1 m   E t + + aE = 0   ρ 2 ρ   x We define the gPC expansions of ρ , m , E and γ in the following form: N N � � ρ N ( x, t, z ) = ρ i ( x, t )Φ i ( z ) , ˆ m N ( x, t, z ) = m i ( x, t )Φ i ( z ) , ˆ i =0 i =0 N N ˆ � � E N ( x, t, z ) = E i ( x, t )Φ i ( z ) , γ N ( z ) = ˆ γ i Φ i ( z ) i =0 i =0 substitute them into the systems and derive the gPC-SG approximation ... 15

We define ... N N 3 − γ N ( z ) ˆ ˆ ˆ � � γ N ( z ) − 1 = ˆ γ i Φ i ( z ) , = γ i Φ i ( z ) , ˆ 2 i =0 i =0 N m 2 � � ˆ � ( x, t, z ) = ψ i ( x, t )Φ i ( z ) , ρ N i =1 N � � γm ˆ ˆ � ( x, t, z ) = ψ i ( x, t )Φ i ( z ) , ρ N i =1 N � � ( γ − 1) m ˆ ˆ ˆ � ( x, t, z ) = ψ i ( x, t )Φ i ( z ) . ρ N i =1 ψ i can be computed by using ρψ = m 2 , namely, For example, ˆ N N ˆ � � ψ k ˆ ρ ℓ S ikℓ = m k ˆ ˆ i = 0 , . . . , N m ℓ S ikℓ , k,ℓ =0 k,ℓ =0 � S ikℓ = Φ i ( z )Φ k ( z )Φ ℓ ( z ) µ ( z ) d z is computed once Ω 16

... after implementing the Galerkin projection we obtain the corresponding three systems for the gPC coefficients i = 0 , . . . , N :  (ˆ ρ i ) t + ( ˆ m i ) x = 0     N    � γ k ( ˆ ˆ ( ˆ m i ) t + ˆ E ℓ ) x S kℓi + ( a ˆ m i ) x = 0 ( I )  k,ℓ =0     ( ˆ E i ) t − ( a ˆ E i ) x = 0    (ˆ ρ i ) t = 0     N   ˆ  ˆ γ k ( ˆ � ( ˆ m i ) t + ˆ ψ ℓ ) x S kℓi − ( a ˆ m i ) x = 0 ( II ) k,ℓ =0      ( ˆ E i ) t = 0    (ˆ ρ i ) t = 0     ( ˆ m i ) t = 0    ( III ) N N ( ˆ ( ˆ ˆ  ˆ ψ k ˆ ˆ ( ˆ � ψ k ˆ � ψ ℓ ) x S kℓi + ( a ˆ E i ) t + E ℓ ) x S kℓi − E i ) x = 0      k,ℓ =0 k,ℓ =0  17