• No results found

A stochastic Galerkin method for the Euler equations with Roe variable transformation

N/A
N/A
Protected

Academic year: 2021

Share "A stochastic Galerkin method for the Euler equations with Roe variable transformation"

Copied!
40
0
0

Loading.... (view fulltext now)

Full text

(1)

A stochastic Galerkin method for the Euler

equations with Roe variable transformation

Per Pettersson, Gianluca Iaccarino and Jan Nordström

Linköping University Post Print

N.B.: When citing this work, cite the original article.

Original Publication:

Per Pettersson, Gianluca Iaccarino and Jan Nordström, A stochastic Galerkin method for the Euler equations with Roe variable transformation, 2013, Journal of Computational Physics.

http://dx.doi.org/10.1016/j.jcp.2013.10.011

Copyright: Elsevier

http://www.elsevier.com/

Postprint available at: Linköping University Electronic Press

(2)

A stochastic Galerkin method for the Euler equations

with Roe variable transformation

Per Petterssona,b,c,∗, Gianluca Iaccarinoa, Jan Nordstr¨omd

aInstitute for Computational and Mathematical Engineering, Stanford University,

Stanford, CA 94305, USA.

bDepartment of Energy Resources Engineering, Stanford University, Stanford, CA 94305,

USA.

cDepartment of Information Technology, Uppsala University, SE-75105 Uppsala, Sweden. dDepartment of Mathematics, Computational Mathematics, Link¨oping University,

SE-58183 Link¨oping, Sweden.

Abstract

The Euler equations subject to uncertainty in the initial and boundary condi-tions are investigated via the stochastic Galerkin approach. We present a new fully intrusive method based on a variable transformation of the continuous equations. Roe variables are employed to get quadratic dependence in the flux function and a well-defined Roe average matrix that can be determined without matrix inversion.

In previous formulations based on generalized polynomial chaos expan-sion of the physical variables, the need to introduce stochastic expanexpan-sions of inverse quantities, or square-roots of stochastic quantities of interest, adds to the number of possible different ways to approximate the original stochastic problem. We present a method where the square roots occur in the choice of variables, resulting in an unambiguous problem formulation.

The Roe formulation saves computational cost compared to the formula-tion based on expansion of conservative variables. Moreover, the Roe formu-lation is more robust and can handle cases of supersonic flow, for which the conservative variable formulation fails to produce a bounded solution. For certain stochastic basis functions, the proposed method can be made more effective and well-conditioned. This leads to increased robustness for both

Corresponding author

Email addresses: massperp@stanford.edu (Per Pettersson ), jops@stanford.edu (Gianluca Iaccarino), jan.nordstrom@liu.se (Jan Nordstr¨om)

(3)

choices of variables. We use a multi-wavelet basis that can be chosen to in-clude a large number of resolution levels to handle more extreme cases (e.g. strong discontinuities) in a robust way. For smooth cases, the order of the polynomial representation can be increased for increased accuracy.

Keywords: Uncertainty quantification, Euler equations, Roe variable

transformation, Stochastic Galerkin method, Multi-wavelets

1. Introduction

Generalized polynomial chaos (gPC) expansions are frequently used to represent uncertain quantities in numerical solutions of partial differential equations (PDEs) with uncertainty in e.g. initial data, boundary data, un-known material parameters and forcing functions due to unun-known geome-try. These quantities of interest are expressed as generalized Fourier series, where orthogonal polynomials (e.g. Legendre or Hermite polynomials) are commonly used as basis functions. Building on results due to Wiener [1], this is the polynomial chaos method introduced by Ghanem and Spanos [2] and generalized in [3].

Spectral convergence of the gPC expansion is observed when the solu-tions are sufficiently regular and continuous [3], but for general non-linear conservation laws - such as in fluid dynamics problems - the convergence is usually less favorable. Spectral expansion representations are still of interest for these problems because of the potential efficiency with respect to brute force sampling methods, but special attention must be devoted to the nu-merical methodology used. For some problems with steep gradients in the stochastic dimensions, gPC expansions fail entirely to capture the solution [4]. Global methods can still give superior overall performance, for instance

Pad´e approximation methods based on rational function approximation [5],

and hierarchical wavelet methods that are global methods with localized sup-port of each resolution level [6, 7].

Intrusive multi-element gPC methods for nonlinear conservation laws have been investigated by Tryoen et. al. in [8], where a reduced-cost Roe solver with entropy corrector was presented, and in [9] with different localized representations of uncertainty in initial functions and problem coefficients. In previous works we investigated well-posedness and time stability for an intrusive formulation of Burgers’ equation with uncertainty in [10], and im-position of uncertain boundary conditions in [11].

(4)

In many nonlinear applications of the stochastic Galerkin method, trun-cation of the gPC expansion leads to non-unique formulations of the systems of equations. For instance, cubic products between stochastic quantities a,

b and c, are represented as products of truncated approximations ˜a, ˜b and

˜

c, but the pseudo-spectral multiplication operator ∗, to be explicitly defined

in a later section, is not associative, i.e. (˜a∗ ˜b) ∗ ˜c 6= ˜a ∗ (˜b ∗ ˜c). Similar problems are investigated in more detail in [12]. It is common practice to introduce these pseudo-spectral approximations since they imply a reduced numerical cost. Examples in the context of polynomial chaos for fluid flow include [8, 9].

The need to introduce gPC expansions of inverse quantities, or square-roots of stochastic quantities of interest, adds to the number of possible different ways to approximate the original stochastic problem. This leads to ambiguity of the problem formulation. We present a method where this ambiguity is avoided. Our formulation relies on a variable transformation where the square-root of the density is computed, a computation that can be performed in a robust way in a small number of operations.

Alternative approaches to gPC methods have also been presented in the literature. Abgrall et. al. [13, 14] developed a semi-intrusive method based on a finite-volume like reconstruction technique of the discretized stochastic space. A deterministic problem is obtained by taking conditional expecta-tions given a stochastic subcell, over which ENO construcexpecta-tions are used to reconstruct the fluxes in the stochastic dimensions. This makes it particu-larly suitable for non-smooth probability distributions.

Po¨ette et. al. [15] used a nonlinear projection method to bound the

os-cillations close to stochastic discontinuities by gPC expansion of the entropy variables obtained from a transformation of the conservative variables. Each time step is complemented by a functional minimization to obtain the en-tropy variables needed to update the solution vector. The method we will present here may appear similar at first sight, but it relies on a different kind of variable transformation and not on kinetic theory considerations. We do not suggest a variable transformation for general conservation laws, but a formulation that specifically targets the solution of the Euler equations with uncertainty in the variables. It is less complicated than a direct gPC expansion of the conservative variables.

In our method, the system of equations is reformulated using Roe vari-ables so that only quadratic terms occur. No fourth-order tensors need to be approximated or calculated, resulting in increased accuracy and reduced

(5)

computational cost. Moreover, there is no need to compute additional chaos expansions for inverse quantities. The Roe variable expansion provides a simple and unambiguous formulation of the Euler equations. For brevity of notation, we will refer to our expansion method as the Roe expansion, and the method based on expansion of the conservative variables as the conservative expansion.

We consider the Sod test case subject to uncertainty in the density, and uncertain diaphragm location, respectively. The uncertainty is represented with a multi-wavelet (MW) expansion in the stochastic dimension, following the framework outlined in [7]. Special cases of the MW basis include the Leg-endre polynomials and the piecewise constant Haar wavelets. The stochastic Galerkin system is obtained by projection of the stochastic Euler equations onto the MW basis functions.

Stochastic hyperbolic problems in general require a large number of stochas-tic basis functions for accurate representation. In parstochas-ticular, this problem becomes severe at large times [16]. A remedy is to use an adaptive stochastic basis that evolves in space and time to save computational cost. In the con-text of stochastic Galerkin methods for hyperbolic problems, Tryoen et. al. introduced an adaptive method where the resolution was determined locally based on numerical estimates of the smoothness of the solution [17]. In this paper, we have chosen to restrict ourselves to a non-adaptive stochastic basis and focus on the numerical solver rather than the stochastic representation. We employ a MUSCL reconstruction in space [18] and a fourth-order Runge-Kutta method for the time integration. A Roe flux based on a stochas-tic Galerkin Roe average matrix will be employed. A Roe average matrix for the numerical Roe flux is derived and proven to fulfill the conditions of Roe [19] under certain circumstances.

In section 2, we present the framework for the stochastic Galerkin for-mulations of the Euler equations in conservative and Roe variables, which are introduced in section 3. In section 4, the numerical flux functions are introduced. We derive and prove properties of the solvers that are necessary to capture essential dynamics of hyperbolic problems. We use a Roe flux suitable for situations where the system eigenvalues can be accurately esti-mated. Numerical results are presented in section 5, where the previously developed techniques are evaluated. Conclusions are drawn in section 6.

(6)

2. Representation of uncertainty

We restrict ourselves to the case of a single random variable in this paper. We believe that the essential features of the numerical method we propose

can be clearly demonstrated in this simplified setting. Let (Ω,A, P) be a

probability space with event space Ω, and probability measure P defined on

the σ-field A of subsets of Ω. Let ξ(ω) be a random variable for ω ∈ Ω.

Consider a generalized chaos basis i(ξ)}∞i=0 spanning the space of second

order (i.e. finite variance) random processes on this probability space. The basis functions are assumed to be orthonormal, i.e. they satisfy

hψiψji = δij, (1)

where the inner product between two functions a(ξ) and b(ξ) is defined by ha(ξ)b(ξ)i =

Z Ω

a(ξ)b(ξ)dP(ξ).

Any second order random field u(x, t, ξ) (i.e. hu2i < ∞) can be expressed as

u(x, t, ξ) = ∞ X

i=0

ui(x, t)ψi(ξ). (2)

Independent of the choice of basis i(ξ)}∞i=0, we can express the mean and variance of u(x, t, ξ) as E(u(x, t, ξ)) = u0(x, t), Var(u(x, t, ξ)) = ∞ X i=1 ui(x, t)2,

respectively. For practical purposes, (2) is truncated to a finite number P terms, and we set

u(x, t, ξ) P X

i=0

ui(x, t)ψi(ξ). (3)

The truncation of the gPC series leading to (3) may result in solutions that are unphysical. An extreme example is when a strictly positive quantity, say density, with uncertainty within a bounded range is represented by a polynomial expansion with infinite range, for instance Hermite polynomials of standard Gaussian variables. The Hermite series expansion converges to

(7)

the true density with bounded range in the limit P → ∞, but for a given order of expansion, say P = 1, the representation ρ = ρ0+ρ1H1(ξ) may result in negative density with non-zero probability since the Hermite polynomial

H1(ξ) = ξ takes arbitrarily large negative values. Similar problems may

be encountered also for polynomial representations with bounded support. Furthermore, polynomial reconstruction of a discontinuity in stochastic space leads to Gibbs oscillations that may yield negative values of an approximation of a solution that is close to zero but strictly positive by definition.

For smooth problems, the gPC expansion is attractive due to its spectral convergence. For non-smooth problems such as nonlinear hyperbolic prob-lems, discontinuities will emerge in finite time in stochastic space and gPC expansions tend to result in Gibbs oscillations. In the worst case, polynomial chaos may fail entirely to capture essential features of the solution [6]. 2.1. Multi-wavelet expansion

An alternative to gPC expansions for non-smooth and oscillatory prob-lems is generalized chaos based on a localization or discretization of the stochastic space [20, 21]. Methods based on stochastic discretization include adaptive stochastic multi-elements [22] and stochastic simplex collocation [23]. The robust properties of discretized stochastic space can also be ob-tained by globally defined wavelets, see [6]. In this paper, we follow the approach of [7] and use piecewise polynomial multi-wavelets (MW), defined

on sub-intervals of [−1, 1]. The construction of a truncated MW basis follows

the algorithm in [24].

Wavelets are defined hierarchically on different resolution levels, repre-senting successively finer features of the solution with increasing resolution. They have non-overlapping support within each resolution level, and in this sense they are localized. Still, the basis is global due to the overlapping support of wavelets belonging to different resolution levels. Piecewise con-stant wavelets, denoted Haar wavelets, do not exhibit spectral convergence, but avoid the Gibbs phenomenon in the proximity of discontinuities in the stochastic dimension.

Starting with the space VNp of polynomials of degree at most Np defined

on the interval [−1, 1], the construction of multi-wavelets aims at finding

a basis of piecewise polynomials for the orthogonal complement of VNp in

the space VNp+1 of polynomials of degree at most Np + 1. Merging the

bases of VNp and that of the orthogonal complement of VNp in VNp+1, we

(8)

finding orthogonal complements in spaces of increasing degree of piecewise polynomials, leads to a basis for L2([−1, 1]).

We first introduce a smooth polynomial basis on [−1, 1]. Let {Lei(ξ)}∞i=0

be the set of Legendre polynomials that are defined on [−1, 1] and orthogonal

with respect to the uniform measure. The normalized Legendre polynomials are defined recursively by

Lei+1(ξ) = √ 2i + 3 √ 2i + 1 i + 1 ξLei(ξ)− i (i + 1)√2i− 1Lei−1(ξ)  , i≥ 1, Le0(ξ) = 1, Le1(ξ) = √ 3ξ. The set {Lei(ξ)}

Np

i=0 is an orthonormal basis for VNp. Double products are

readily computed from (1), and higher-order products are precomputed using numerical integration.

Following the algorithm by Alpert [24] (see Appendix A), we construct a

set of mother wavelets W

i (ξ)}

Np

i=0 defined on the domain ξ∈ [−1, 1], where

ψiW(ξ) =    pi(ξ) −1 ≤ ξ < 0 (−1)Np+i+1p i(ξ) 0≤ ξ < 1 0 otherwise, (4)

where pi(ξ) is an ith order polynomial. By construction, the set of wavelets

{ψW

i (ξ)}

Np

i=0 are orthogonal to all polynomials of order at most Np, hence

the wavelets are orthogonal to the set {Lei(ξ)}

Np

i=0 of Legendre polynomials

of order at most Np. Based on translations and dilations of (4), we get the

wavelet family

ψi,j,kW (ξ) = 2j/2ψWi (2jξ− k), i = 0, ..., Np, j = 0, 1, ..., k = 0, ..., 2j−1.

Let ψm+1(ξ) for m = 0, . . . , Np be the set of Legendre polynomials up to

order Np, and concatenate the indices i, j, k into m = (Np+1)(2j+k−1)+i+1

so that ψm(ξ) ≡ ψi,j,kW (ξ) for m > Np + 1. With the MW basis {ψm(ξ)}∞m=0 we can represent any random variable u(x, t, ξ) with finite variance as

u(x, t, ξ) = ∞ X m=1

um(x, t)ψm(ξ),

which is of the form (2), but with index starting from 1 in order to simplify the notation henceforth. In the computations, we truncate the MW series

(9)

both in terms of the piecewise polynomial order Np and the resolution level

Nr. With the index j = 0, ..., Nr, we retain P = (Np+ 1)2Nr terms of the

MW expansion.

The truncated MW basis is characterized by the piecewise polynomial

order Np and the number of resolution levels Nr, illustrated in Figure 1 for

Np = 2 and Nr = 3. As special cases of the MW basis, we obtain the

Legendre polynomial basis for Nr = 0 (i = j = 0), and the Haar wavelet

basis of piecewise constant functions for Np = 0.

−1 −0.5 0 0.5 1 −2 −1 0 1 2 Resolution level 0 −1 −0.5 0 0.5 1 −4 −2 0 2 4 Resolution level 1 −1 −0.5 0 0.5 1 −5 0 5 Resolution level 2 −1 −0.5 0 0.5 1 −5 0 5 Resolution level 3

Figure 1: Multi-wavelets for Np= 2, Nr= 3. Resolution level 0 consists of the first Np+ 1

Legendre polynomials. Resolution level j > 0 contains 2j−1 wavelets each. Each basis

function is a piecewise polynomial of order Np.

3. Euler equations with input uncertainty

Consider the 1D Euler equations, in non-dimensional form given by

ut+ f (u)x = 0, 0≤ x ≤ 1, t > 0, (5)

where the solution and flux vector are given by

u =   ρ ρv E   , f =   ρv ρv2+ p (E + p)v   ,

(10)

where ρ is density, v velocity, E total energy and pressure p. A perfect gas equation of state is assumed, and energy and pressure are related by

E = p

γ− 1 +

1 2ρv

2,

where γ is the ratio of the specific heats. For the numerical method, we need the flux Jacobian, given by

∂f ∂u =   0 1 0 1 2(γ− 3)v 2 (3− γ)v γ− 1 1 2(γ− 1)v 3− vH H − (γ − 1)v2 γv   , with the total enthalpy H = (E + p)/ρ.

We scale the physical variables to get the dimensionless variables ρ = ρ0/ρ0ref, E = E0/(γp0ref), p = p0/(γpref0 ) and v = v0/a0ref where a0 = (γp0/ρ0)1/2 and the subscript ref denotes a reference state.

3.1. Formulation in Roe variables Roe [19] introduced the variables

w =   w1 w2 w3   =   ρ1/2 ρ1/2v ρ1/2H   . The flux and the conservative variables are given by

ˆ f (w) =   w1w2 γ−1 γ w1w3+ γ+1 2γ w 2 2 w2w3   , u = ˆg(w) =   w2 1 w1w2 w1w3 γ + γ−1 2γ w 2 2   . Then ˆ g(w)t+ ˆfx(w) = 0 (6)

is equivalent to (5). The flux Jacobian in the Roe variables is given by ∂ ˆf ∂w =   w2 w1 0 γ−1 γ w3 γ+1 γ w2 γ−1 γ w1 0 w3 w2   .

(11)

3.2. Stochastic Galerkin formulation of the Euler equations

Define the pseudo-spectral product u∗ v of order P = P (Np, Nr) by

(u∗ v)k= P X i=1 P X j=1 uivjhψiψjψki, k = 1, . . . , P, where hψiψjψki = Z Ω ψi(ξ)ψj(ξ)ψk(ξ)dP.

Alternatively, using matrix notation, we can write the spectral product as

u∗ v = A(u)v, where [A(u)]jk = P X i=1 uihψiψjψki, j, k = 1, . . . , P. (7)

We will need the pseudo-spectral inverse q−∗, defined as the solution of q

q−∗ = 1, and the pseudo-spectral square root, defined as the solution q∗/2 of

q∗/2∗ q∗/2 = q, where the spectral expansion of the quantity of interest q is

assumed to be known. For more details, see [25].

Let uP denote the vector of coefficients of the MW expansion of u of

order P = P (Np, Nr). P may take the same value for two distinct pairs

of (Np, Nr) but this ambiguity in notation will not matter in the derivation

of the numerical method so for brevity we use only P in the superscripts. The Euler equations represented by the conservative formulation (5) can be written as an augmented system, after stochastic Galerkin projection,

uPt + fP(uP)x = 0, (8) where uP =   uP 1 uP 2 uP3   =   [(u1)0, ..., (u1)P]T [(u2)0, ..., (u2)P] T [(u3)0, ..., (u3)P]T   , fP(uP) =   uP 2 (uP 1) −∗ ∗ uP 2 ∗ uP2 + pP (uP3 + pP)∗ uP2 ∗ (uP1)−∗   . with pP = (γ− 1)(uP 3 − (uP1) −∗∗ uP

2 ∗ uP2/2). The cubic products of (8) are approximated by the application of two third-order tensors, instead of one

fourth-order tensor. That is, we replace (a∗b∗c)l=

P

ijkhψiψjψkψlaibjcki by

(12)

error in addition to the error from the truncation of the gPC series to a finite number of terms. The effect of the error introduced by the approximation of higher order tensors with successive application of third order tensors was studied in [25], where it was found that the error is negligible if sufficiently high order gPC expansions are used. We use this approximation with the conservative variables to make a fair comparison of the computational cost with the method we propose based on Roe variables.

For the Roe variable formulation, the stochastic Galerkin projection of (6) gives the system

ˆ gP(wP)t+ ˆfP(wP)x = 0, (9) where ˆ gP(wP) =   wP 1 ∗ w1P wP1 ∗ w2P wP 1∗wP3 γ + γ−1 2γ w P 2 ∗ w2P   , ˆfP(wP) =   wP 1 ∗ wP2 γ−1 γ w P 1 ∗ w3P + γ+1 2γ w P 2 ∗ wP2 wP 2 ∗ wP3   . The flux Jacobian for the stochastic Galerkin system in the Roe variables is given by ∂ ˆfP ∂wP =   A(wP 2) A(wP1) 0P ×P γ−1 γ A(w P 3) γ+1 γ A(w P 2) γ−1 γ A(w P 1) 0P ×P A(wP3) A(w2P)   . (10)

As P → ∞, the formulations (8) and (9), as well as any other consistent

formulation, are equivalent. However, P is assumed to be small (< 20), and truncation and conditioning of the system matrices will play an important role for the accuracy of the solution.

We assume that γ is a deterministic constant in the formulation of the numerical schemes. Although it would imply additional pseudo-spectral mul-tiplications, accounting for a random γ is a straightforward extension of the

presented framework. This amounts to forming A(γ−1), which can be

pre-computed and stored for use in the updates of the numerical fluxes. 4. Numerical method

We use the MUSCL (Monotone Upstream-centered Schemes for Conser-vation Laws) scheme introduced in [18]. For clarity of comparison of the numerical results, the MUSCL scheme is used both for the conservative vari-able formulation and for the Roe varivari-able formulation.

(13)

4.1. Expansion of conservative variables

Let m be the number of spatial points and ∆x = 1/(m− 1) and let UP

be the spatial discretization of uP. The semi-discretized form of (8) is given

by dUP j dt + FP j+1/2− F P j−1/2 ∆x = 0, j = 1, ..., m, (11)

where Fj+1/2P denotes the numerical flux function evaluated at the interface

between cells j and j + 1.

For the MUSCL scheme with slope limited states UL and UR, we take

the numerical flux Fj+P 1 2 = 1 2  fP(Uj+L 1 2 ) + fP(Uj+R1 2 )+1 2|( ˜J P c )j+12|  Uj+L 1 2 − U R j+12  , (12)

where the Roe average ˜JP

c is the pseudo-spectral generalization of the

stan-dard Roe average of the deterministic Euler equations, i.e.

˜ JcP(v, H) =   0P ×P IP ×P 0P ×P 1 2(γ− 3)A(v) 2 (3− γ)A(v) − 1)I P ×P 1 2(γ− 1)A(v)

3− A(v)A(H) A(H) − (γ − 1)A(v)2 γA(v)

  where v = (ρ−∗/2L + ρ−∗/2R )∗ (ρ∗/2L ∗ vL+ ρ ∗/2 R ∗ vR), and H = (ρ∗/2L ∗ HL+ ρ ∗/2 R ∗ HR)∗ (ρ −∗/2 L + ρ −∗/2 R ).

The computation of v and H require the spectral square root ρ∗/2 and its

inverse, that are computed solving a nonlinear and a linear system, respec-tively.

Further details about the formulation of the Roe average matrix are given in [9]. The scheme is a direct generalization of the deterministic MUSCL scheme. Flux limiters are applied componentwise to all MW coefficients in sharp regions. For a more detailed description of the MUSCL scheme and flux limiters, see e.g. [26], and for application to the stochastic Burgers’ equation [27].

(14)

4.2. Expansion of Roe’s variables

Let WP denote the spatial discretization of wP. The semi-discretized

form of (9) is given by ∂ ˆgP j ∂t + ˆ FP j+1/2− ˆFj−1/2P ∆x = 0, j = 1, ..., m,

with the numerical flux function ˆ Fj+1 2 = 1 2  ˆ fP(Wj+L 1 2 ) + ˆfP(Wj+R1 2 )  +1 2| ˜J P j+12|  Wj+L 1 2 − W R j+12  , (13)

where ˜JP = ˜JP(WP) is the Roe matrix for the stochastic Galerkin formula-tion of the Euler equaformula-tions in Roe’s variables, to be derived below.

Each time step provides the update of the solution vector ˆgP

j , j = 1, ..., m,

from which we can solve for WP to be used in the update of the numerical

flux. This involves solving the nonlinear systems

A(W1,jP )W1,jP = ˆgP1,j, j = 1, ..., m, (14)

for W1,jP , and then using W1,jP to solve the linear P × P -systems A(W1,jP )W2,jP = ˆgP2,j, j = 1, ..., m, for WP 2,j, and A(W1,jP )W3,jP = γ ˆg3,jP γ− 1 2 A(W P 2,j)W2,jP , j = 1, ..., m, for WP 3,j.

The system (14) is solved iteratively with a trust-region-dogleg

algo-rithm1. Starting with the value of the previous time-step as initial guess,

few iterations are required (typically 2-3). (The same method is used to solve for spectral square roots in the conservative variable formulation.)

(15)

4.3. Stochastic Galerkin Roe average matrix for Roe variables

The Roe average matrix ˜JP is given as a function of the Roe variables

w = (w1, w2, w3)T, where each wi is a vector of generalized chaos coefficients. It is designed to satisfy the following properties:

(i) ˜JP(wL, wR) ∂ ˆfP ∂w w=w0 as w L, wR→ w0. (ii) ˜JP(wL, wR)× (wL− wR) = ˆfP(wL)− ˆfP(wR),∀wL, wR

(iii) ˜JP is diagonalizable with real eigenvalues and linearly independent

eigenvectors.

In the standard approach introduced by Roe and commonly used for deter-ministic calculations, the conservative variables are mapped to the w vari-ables, which are then averaged.

In the deterministic case, we have ˆ fL− ˆfR = ˜J (wL, wR)× (wL− wR), (15) where ˜ J (wL, wR) =   w2 w1 0 γ−1 γ w3 γ+1 γ w2 γ−1 γ w1 0 w3 w2   .

Overbars denote arithmetic averages of assumed left and right values of a variable, i.e.

wj =

wL

j + wjR

2 , j = 1, 2, 3.

It is a straightforward extension of the analysis by Roe in [19] to show properties (i) and (ii) for the Roe variables, without mapping to the con-servative variables. To prove (iii) we note that there exists an eigenvalue decomposition ˜ J = V DV−1, (16) where V =    w1 w3 w1 w3 − w1 w3 w2−√w22+8w1w3γ(γ−1) 2γw3 w2+√w22+8w1w3γ(γ−1) 2γw3 0 1 1 1    , (17)

(16)

D =     w2(1+2γ)− √ 8w1w3γ(γ−1)+w22 2γ 0 0 0 w2(1+2γ)+ √ 8w1w3γ(γ−1)+w22 2γ 0 0 0 w2    . (18)

The eigenvalues of ˜J are real and distinct, so property (iii) is also satisfied. Now consider the stochastic Galerkin formulation, i.e. assume that the

wi’s are vectors of generalized chaos coefficients. The stochastic Galerkin

Roe average matrix ˜JP for the Roe variables formulation is a generalization

of the mapping (15), i.e. of the matrix ˜J . We define

˜ JP(wL, wR) = ˜JP(w) =   A(w2) A(w1) 0P ×P γ−1 γ A(w3) γ+1 γ A(w2) γ−1 γ A(w1) 0P ×P A(w3) A(w2)   , (19)

where the submatrix A(wj) is given by (7) and w = (wL+ wR)/2.

Proposition 1. Property (i) is satisfied by (19).

Proof. With wL = wR = w0, ˜JP(wL, wR) = ˜JP(w0, w0) = ∂ ˆfP ∂wP w=w0 by (10).

Proposition 2. Property (ii) is satisfied by (19). Proof. ˜ JP(wL, wR)× (wL− wR) = 1 2  ˜ JP(wL) + ˜JP(wR)(wL− wR) = = 1 2 ˜ JP(wL)wL 1 2 ˜ JP(wR)wR= ˆfP(wL)− ˆfP(wR), (20) where the last equality follows from the fact that the stochastic Galerkin generalizations of the Euler equations are homogeneous of degree 1.

To prove (iii), we will need the following lemma.

Lemma 1. Let A(wj) (j = 1, 2, 3) be defined by (7) and A(wj) = QΛjQT be

an eigenvalue decomposition with constant eigenvector matrix Q and assume

that Λ1 and Λ3 are non-singular. Then the stochastic Galerkin Roe average

matrix ˜JP has an eigenvalue decomposition ˜JP = X ˜ΛPX−1 with a complete set of eigenvectors.

(17)

Proof. We will use the Kronecker product⊗, defined for two matrices B (of size m× n) and C by B⊗ C =    b11C . . . b1nC .. . . .. ... bm1C . . . bmnC    .

The eigenvalue decompositions of each P × P matrix block of (19) have the

same eigenvector matrix Q, hence we can write ˜ JP = (I3×3⊗ Q) ˆJ (I3×3⊗ QT) (21) where ˆ J =   Λ2 Λ1 0P ×P γ−1 γ Λ3 γ+1 γ Λ2 γ−1 γ Λ1 0P ×P Λ3 Λ2   .

By assumption, I3×3⊗ Q is non-singular, and it remains to show that ˆJ has

distinct eigenvectors. Let S = diag(Λ1Λ−13 ,

p

− 1)/γΛ1/21 Λ−1/23 , IP ×P).

By assumption, Λ1 and Λ3 are invertible, so S and S−1 exist. We have

JS ≡ S−1J S =ˆ       Λ2 h γ−1 γ Λ1Λ3 i1/2 0P ×P h γ−1 γ Λ1Λ3 i1/2 γ+1 γ Λ2 h γ−1 γ Λ1Λ3 i1/2 0P ×P h γ−1 γ Λ1Λ3 i1/2 Λ2      . (22)

Clearly, JS is symmetric and has the same eigenvalues as ˆJ and ˜JP. Hence,

JS has an eigenvalue decomposition JS = Y ˜ΛPYT. Then,

ˆ

J = SY ˜ΛPYTS−1 = SY ˜ΛP(SY )−1. (23)

Combining (21) and (23), we get ˜

JP = [(I3×3⊗ Q)SY ]˜ΛP[(I3×3⊗ Q)SY ]−1.

Setting X = (I3×3 ⊗ Q)SY , we get the eigenvalue decomposition ˜JP =

X ˜ΛPX−1. By assumption, S and Y are non-singular, and we have that

det (X) = det ((I3×3⊗ Q)SY ) 6= 0,

which proves that X is non-singular, and thus ˜JP has a complete set of

(18)

Proposition 3. Property (iii) is satisfied by (19).

Proof. Lemma 1 shows that since the eigenvalue matrix ˜ΛP is also the

eigen-value matrix of the symmetric matrix JS defined in (22), the eigenvalues are

all real. Lemma 1 also shows that the eigenvectors are distinct.

The conditions in Lemma 1 are true for certain basis functions assuming moderate stochastic variation, but it can not be guaranteed for every case, and it certainly does not hold for pathological cases with e.g. negative den-sity. The requirement of non-singularity of Λ1, Λ3 is not very restrictive since it amounts to excluding unphysical behavior, for instance naturally positive quantities taking negative values with non-zero probability. The assump-tion of constant eigenvectors of the matrix A holds for Haar wavelets (i.e.

multi-wavelets with Np = 0), for all orders P = 2Nr, with Nr ∈ N. See

appendix B for a proof sketch. Expressions for the first constant eigenvalue decompositions are included in appendix C for Haar wavelets and piecewise linear multi-wavelets. The eigenvectors of A for P = 1, 2, 4, 8 are shown to be constant, but we do not give a proof that this is true for piecewise linear multi-wavelets of any order P .

Remark 1. The Roe variable scheme has been outlined under the implicit assumption of uncertainty manifest in the variables, e.g. initial and bound-ary condition uncertainty. However, situations such as uncertainty in the adiabatic coefficient γ may be treated in a similar way, although it would re-sult in additional pseudo-spectral products. Pseudo-spectral approximations

of (γ− 1)/γ and (γ + 1)/γ could then be precomputed to sufficient accuracy.

Remark 2. For both the conservative variable formulation and the Roe vari-able formulation, we need to find the eigenvalue decomposition of ˜JP

c ( or ˜JP) at each time step and each spatial point. For the cases of piecewise constant or piecewise linear MW, we can do this analytically and thus at low com-putational cost. For higher-order polynomial MW, we may rely on iterative

methods for the eigenvalue decomposition of these 3P × 3P subsystems. To

this end, one may use e.g. the approximate low-order polynomial method, that was introduced and successfully applied in [9] for very similar problems. 5. Numerical results

We use the method of manufactured solutions to verify the second order convergence in space of a smooth problem using the MUSCL scheme with Roe

(19)

variables. We then introduce two test cases for the non-smooth problem; case 1 with an initial function that can be exactly represented by two Legendre polynomials, and case 2 with slow initial decay of the MW coefficients in both

Np and Nr. The errors in computed quantities of interest (here variances)

as functions of the order of MW are investigated. Qualitative results are then presented to indicate the behavior we can expect for the convergence of two special cases of MW, namely the Legendre polynomials and Haar wavelet basis, respectively. Robustness with respect to more extreme cases (density close to zero leading to high Mach number) is demonstrated for the Roe variable formulation for a supersonic case where the conservative variable method breaks down. Finally, we perform a comparative study of the computational time for the formulation in conservative variables and the formulation in Roe variables.

5.1. Spatial convergence

The MUSCL scheme with appropriate flux limiters is second order ac-curate for smooth solutions. Since the Euler solution in general becomes discontinuous in finite time, the method of manufactured solutions [29, 30] is used to solve the Euler equations with source terms for a known smooth solu-tion. The smooth solution is inserted into the Euler equations (5) and results in a non-zero right-hand side that is used as a source function. In order to test the capabilities of the method, we choose a solution that varies in space, time and in the stochastic dimension and with time-dependent boundary conditions. It is designed to resemble a physical solution with non-negative density and pressure. The solution is given by

  ρ v p   =   ρ0+ ρ1tanh(s(x0− x + t + σξ)) tanh(s(x0+ v0− x + t + σξ)) + tanh(−s(x0− v0− x + t + σξ)) p0+ p1tanh(s(x0− x + t + σξ))   .

The parameters are set to ρ0 = p0 = 0.75, ρ1 = p1 = x0 = 0.25, v0 = 0.05,

s = 10, σ = 0.1 and ξ ∈ U[−1, 1].

(20)

and the discrete `2 norm, uP − u 2,2 ≡ uP − u `2,L2(Ω,P) = = ∆x m X i=1 uP(x i, t, ξ)− u(xi, t, ξ) 2 L2(Ω,P) !1/2 = = ∆x m X i=1 Z Ω (uP(xi, t, ξ)− u(xi, t, ξ))2dP(ξ) !1/2 = ≈ ∆x m X i=1 q X j=1 (uP(xi, t, ξq(j))− u(xi, t, ξ(j)q )) 2w(j) q !1/2 , (24)

where a q-point quadrature rule with points q(j)}qj=1 and weigths {w (j) q }qj=1

was used in the last line to approximate the integral in ξ. The

Gauss-Legendre quadrature is used here since the solution is smooth in the stochas-tic dimension.

Figure 2 depicts the spatial convergence in the k.k2,2 norm of the error

in density, velocity and energy. An order (Np, Nr) = (10, 0) basis is used to represent the uncertainty. The solution dynamics is initially concentrated in the left part of the spatial domain. By the time of t = 0.4, it has moved to the right and has begun to exit the spatial domain, so the time snapshots of Figure 2 summarizes the temporal history of the spatial error decay. The theoretical optimal convergence rate for the MUSCL scheme with the van Leer flux limiter is obtained for all times and all quantities.

(21)

21 41 81 161 10−4 10−3 10−2 10−1 m ° °ρP− ρ°° 2 , 2 ° °vP− v°° 2 , 2 ° °EP− E°° 2 , 2 or d . 2 d e c . (a) t = 0.05. 21 41 81 161 10−4 10−3 10−2 10−1 m ° °ρP− ρ°° 2 , 2 ° °vP− v°° 2 , 2 ° °EP− E°° 2 , 2 or d . 2 d e c . (b) t = 0.1. 21 41 81 161 10−4 10−3 10−2 10−1 m ° °ρP− ρ°° 2 , 2 ° °vP− v°° 2 , 2 ° °EP− E°° 2 , 2 or d . 2 d e c . (c) t = 0.2. 21 41 81 161 10−4 10−3 10−2 10−1 m ° °ρP− ρ°° 2 , 2 ° °vP− v°° 2 , 2 ° °EP− E°° 2 , 2 or d . 2 d e c . (d) t = 0.4.

Figure 2: Convergence in space using the method of manufactured solutions, Np = 10,

Nr= 0 (Legendre polynomials).

5.2. Initial conditions and discontinuous solutions

We consider (5) with two different initial functions on the domain [0, 1]. Since the analytical solution of Sod’s test case is known for any fixed value of the input parameters, the exact stochastic solution can be formulated as a function of the stochastic input ξ. Exact statistics can be computed by numerical integration over ξ. As case 1, assume that the density is subject to uncertainty, and all other quantities are deterministic at t = 0. The initial

(22)

condition for (5) is given by u(x, t = 0, ξ) =



uL= (1 + σξ, 0, 2.5/γ)T x < 0.5

uR= (0.125(1 + σξ), 0, 0.25/γ)T x > 0.5

where we assume ξ ∈ U[−1, 1], γ = 1.4 and the scaling parameter σ = 0.5.

This is a simple initial condition in the sense that the first two Legendre polynomials are sufficient to represent the initial function exactly. As case no 2, we consider (5) subject to uncertainty in the initial shock location. Let

u(x, t = 0, ξ) = 

uL = (1, 0, 2.5/γ)T x < 0.5 + ση

uR = (0.125, 0, 0.25/γ)T x > 0.5 + ση

where we assume γ = 1.4 and the scaling parameter σ = 0.05. Here, η takes a triangular distribution, which we parameterize as a nonlinear function in ξ ∈ U[−1, 1], i.e.

η(ξ) = (−1 +pξ + 1)1{−1≤ξ≤0}(ξ) + (1−

p

1− ξ)1{0<ξ≤1}(ξ),

where the indicator function 1{A} of a set A is defined by 1{A}(ξ) = 1 if

ξ ∈ A and zero otherwise. For case 2, exact representation of the initial

function requires an infinite number of expansion terms in the MW basis. Figure 3 depicts the shock tube setup for the two cases, with dashed lines denoting uncertain parameters. We will also investigate another version of case 2, where the right state density is significantly reduced to obtain a strong shock. 0 x0 1 ρL ρR 2σ ? 6 0.25σ ? 6 2σ  -0 x0 1 ρL ρR

Figure 3: Schematic representation of the initial setup for case 1 (left) and case 2 (right).

5.3. Initial conditions and resolution requirements

For case 2, it should be noted that although the initial shock position can be exactly described by the first two terms of the Legendre polynomial chaos

(23)

expansion, this is not the case for the initial state variables. In fact, for the the polynomial chaos expansions of the density, momentum and energy, the error decay only slowly with the number of expansion terms. Thus, unless a reasonably large number of expansion terms are retained, the stochastic Galerkin solution of case 2 will not be accurate even for small times.

The Legendre coefficients at small times display an oscillating behavior that becomes sharper with the order of the coefficients. The wavelet coeffi-cients exhibit peaks that get sharper with the resolution level, and require a fine mesh. Figure 4 shows the initial Legendre coefficients and the initial Haar wavelets for case 2. The numerical method has a tendency to smear the chaos coefficients, resulting in under-predicting the variance. The increasing cost of using a larger number of basis functions is further increased by the need for a finer mesh to resolve the solution modes.

Figure 5 shows the temporal evolution of the mean and variance of the density of case 2 as a function of space on a fine mesh of 500 spatial points and 32 piecewise linear multi-wavelets. The mean and the variance are both reasonably well captured for this case. Figure 6 depicts case 2 for a similar setup, but with 32 Haar wavelets. The mean is well captured whereas the variance is not fully captured. The three variance peaks correspond to the rarefaction wave, contact discontinuity and shock, respectively. As time pro-gresses, the variance peaks will propagate out of the computational domain.

0.3 0.4 0.5 0.6 0.7 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 (w 1)0 (w 1)1 (w 1)2 (w 1)3 (w 1)4 (w 1)5 (w 1)6 (w 1)7

(a) Legendre polynomials.

0.3 0.4 0.5 0.6 0.7 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 (w 1)0 (w 1)1 (w 1)2 (w 1)3 (w 1)4 (w 1)5 (w 1)6 (w 1)7 (b) Haar wavelets. Figure 4: Initial w1modes for case 2, first 8 basis functions.

(24)

0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 MW, t=0.1 Ref, t=0.1 MW, t=0.3 Ref, t=0.3

(a) Mean density.

0 0.2 0.4 0.6 0.8 1 0 0.02 0.04 0.06 0.08 0.1 0.12 MW, t=0.1 Ref, t=0.1 MW, t=0.3 Ref, t=0.3 (b) Variance of density.

Figure 5: Temporal evolution of the mean and variance of the density for case 1, using Roe variables, 500 spatial points and 32 piecewise linear multi-wavelets.

0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 MW, t=0.1 Ref, t=0.1 MW, t=0.3 Ref, t=0.3

(a) Mean density.

0 0.2 0.4 0.6 0.8 1 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 MW, t=0.1 Ref, t=0.1 MW, t=0.3 Ref, t=0.3 (b) Variance of density.

Figure 6: Temporal evolution of the mean and variance of the density for case 2, using Roe variables, 500 spatial points and 32 Haar wavelets.

5.4. Convergence of multi-wavelet expansions

For moderate simulation times, the numerical solution on a sufficiently fine spatial mesh converges as the order of MW expansion increases by

(25)

the decay in the error of the variance of velocity and energy as a function of

Np and Nr. For well-behaved cases like these, one may freely choose between

increasing Np and Nr, in order to increase the accuracy of the solution of the quantity of interest. 0 1 2 3 1 2 3 0 0.01 0.02 0.03 0.04 N r N p

(a) Case 1, V ar(vP)

− V ar(v) 2. 0 1 2 3 1 2 3 0 0.002 0.004 0.006 0.008 0.01 N r N p (b) Case 1, V ar(EP) − V ar(E) 2. 0 1 2 3 1 2 3 0 0.01 0.02 0.03 0.04 N r N p

(c) Case 2, V ar(vP)− V ar(v) 2. 0 1 2 3 1 2 3 0 0.01 0.02 0.03 0.04 N r N p

(d) Case 2, V ar(EP)− V ar(E) 2.

Figure 7: Decay in variance of velocity and energy as a function of the order of expansion, polynomial order Npand resolution level Nr. Case 1, t = 0.05, 280 spatial points restricted

(26)

For longer simulation times or more extreme cases, e.g. supersonic flow,

high-order polynomial representation (increasing Np) may not lead to

in-creased accuracy, but instead breakdown of the numerical method. Next, we study the qualitative properties of the MW representation of case 1 and

case 2 for two extreme cases: Legendre polynomials (Nr = 0) and piecewise

constant Haar wavelets (Np = 0).

Figure 8 shows the density surface in the x− ξ-plane of case 1 and case 2

at t = 0.15 based on exact solution evaluations, and computed with Legen-dre polynomials and Haar wavelets. The computed solution with LegenLegen-dre polynomial reconstruction captures essential features of the exact solution, but the use of global polynomials cause oscillations downstream of the shock. With Haar wavelets, there are no oscillations downstream, as in the Leg-endre polynomials case. However, the 8 ’plateaus’ seen in figure 8 (e) corre-sponds to the 8 basis functions. When the order of wavelet chaos expansion increases, the number of plateaus increases, and the solution converges to the exact solution.

From Figure 8 it is clear that the effect of the choice of multi-wavelet basis to some extent depends on the problem at hand. The Haar wavelets yield numerical solutions that are free of oscillations but converge only slowly. Oscillations around discontinuities in stochastic space should be expected when a polynomial basis is used and may lead to severe problems when variables attain unphysical values, e.g. when the oscillations downstream of the shock leads to negative density. Thus, more robust multi-wavelets are required for problems with stronger shocks, as we demonstrate below. 5.5. Robustness

Complex supersonic test cases have already been successfully treated with a stochastic Galerkin method based on the conservative formulation, see for instance [8]. In general, the stochastic Galerkin method applied to the Roe variables gives a more robust method than the conservative variables formu-lation. The conservative formulation is more prone to ill-conditioning of the pseudo-spectral operations in the computation of the numerical flux. How-ever, for the cases where the matrix A in (7) has an eigenvalue decomposition with constant eigenvectors, the pseudo-spectral systems simplify to a series of scalar operations, thus avoiding ill-conditioned systems. An example is given in the Section 5.6.

Figure 9 shows the relative errors of the solution in the 2, 2 norm for

(27)

0 0.5 1 −1 0 1 0 0.5 1 1.5 ξ x

(a) Exact solution, case 1.

0 0.5 1 −1 0 1 0 0.2 0.4 0.6 0.8 1 ξ x

(b) Exact solution, case 2.

0 0.5 1 −1 0 1 0 0.5 1 1.5 ξ x (c) 8 Legendre polynomials (Np, Nr) = (8, 0), case 1. 0 0.5 1 −1 0 1 0 0.2 0.4 0.6 0.8 1 ξ x (d) 8 Legendre polynomials (Np, Nr) = (8, 0), case 2. 0 0.5 1 −1 0 1 0 0.5 1 1.5 ξ x

(e) Haar wavelets (Np, Nr) = (0, 3), case 1. 0 0.5 1 −1 0 1 0 0.2 0.4 0.6 0.8 1 ξ x

(f) Haar wavelets (Np, Nr) = (0, 3), case 2.

(28)

state densities, ρR = 2−k, k = 3, .., 8 for 8 basis wavelets. This corresponds to Mach numbers up to M a = 2.0. Figure 9 also includes the relative error of the Mach number to verify that the cases solved for were reasonable close to the supersonic range they model. For this problem, the conservative vari-able formulation was unstvari-able due to ill-conditioning of the pseudo-spectral

operations except for the original subsonic case 2 (ρR = 0.125). Note that

this numerical breakdown should not be confused with time instability - us-ing analytical decomposition of the eigenvectors of A defined in (7), we can handle supersonic cases also with the conservative variables. We have not observed any significant variation in the stability properties depending on the order P of the stochastic basis when using constant eigenvector decom-positions. Thus, although no clear explanation is known to us, it seems that the Roe variable formulation is more suitable for problems where robustness is an issue, unless the eigenvectors of A are constant.

Legendre polynomials are not suitable for this problem. As seen in Figure 8 (c) and (d), the solution is oscillatory in the right state close to the shock. If the right state density is small, as in this supersonic case, such oscillations cause the density to be very close to zero, or even negative. This leads to an unphysical solution and breakdown of the numerical method.

−3 −4 −5 −6 −7 −8 0 0.02 0.04 0.06 0.08 0.1 0.12 log2( ρR) kρP− ρk 2 , 2 kρk2 , 2 kvP − vk2 , 2 kv k2 , 2 kEP − Ek2 , 2 kE k2 , 2 kM aP− M a k2 , 2 kM ak2 , 2

Figure 9: Relative error in density, velocity, energy and Mach number at t = 0.15 for different shock strengths. m = 300 spatial points, 8 Haar wavelets (Np= 0, Nr= 3).

(29)

5.6. Computational cost

For stochastic basis functions that admit an eigenvalue decomposition of the matrix A in (7) with constant eigenvectors, the computational cost is greatly reduced compared to the general case of non-constant eigenvectors. The Roe average matrices are computed by a series of matrix-vector mul-tiplications only, both for the Roe variables and the conservative variables. The nonlinear pseudo-spectral operations are also simplified. For instance, the computation of the pseudo-spectral inverse used with the conservative variables can be computed by a series of scalar inverses and matrix-vector multiplications. Let Q be the matrix of constant eigenvectors of A(.).

Start-ing from the gPC expansion ρP of the density, put ρ−∗ = 1/P Qρ−∗

EV, where

the vector ρ−∗EV is defined by (ρ−∗EV)j = 1/( √

P QTρ)

j for j = 1, . . . , P . To see

that this holds, note that ΛPρ =√P diag(QTρP). (Superscript P denotes an

index, not a power.) Then

A(ρ)ρ−∗ = QΛPρQTQρ−∗EV = QΛPρρ−∗EV = Q1P = e1,

so ρ−∗ has the desired properties of the pseudo-spectral inverse.

For two stochastic Galerkin systems of order P = (Np + 1)2Nr and

P0 = (Np0 + 1)2Nr0 where P = P0 but N

p 6= Np0, Nr0 6= Nr, the size of the problem and the computational cost is the same. Although the different bases could possibly result in properties that make them very different in e.g. the number of iterations required to solve the nonlinear matrix prob-lems, no such tendency was observed. The numerical experiments yield very similar computational costs for the cases tested.

In order to compare the computational cost of the Roe variable expansion with that of the conservative expansion, a similar experimental setup is used for both methods. Sufficiently small test cases are run in order not to exceed the cache limit which would slow down the simulation time for fine meshes and bias the result. We used test case 1 for short simulation times. Results are shown both for the numerical methods where we use the knowledge of the constant eigenvectors of the eigenvalue decomposition of A, and for the methods designed for the more general case of varying eigenvectors where we have to rely on methods for nonlinear systems.

In the experiments, the same time step has been used for the different variable expansions although a larger time step could be used for the Roe variables. Table 1 displays the relative simulation time of the two different

variable expansions for increasing number of Haar wavelets (P = 2Nr, N

(30)

0). One time unit is defined as the time for the numerical simulation of a single deterministic problem using the same numerical method with similar input conditions, discretization and time step. In the general setup, the higher computational cost for the conservative variable formulation is due to the need to compute inverse quantities and cubic spectral products. The Roe variable formulation only requires the solution of the nonlinear system for the square root of the density and quadratic flux function evaluations. The relative benefit of the Roe variable expansion decreases with the order of wavelet expansion. This is due to the increasing cost of forming spectral products that dominates the total cost for high-order expansions. Note that the difference in computational cost is too large to be due only to the fact that additional pseudo-spectral operations are required for the conservative variables. The difference is attributed to the fact that the ill-conditioned pseudo-spectral operations may need a large number of iterations.

Using the constant eigenvectors of the eigenvalue decomposition of A, we see in Table 1 that the two formulations are essentially equivalent since they both reduce to a comparable number of matrix-vector products instead of solution of nonlinear and possibly ill-conditioned systems. Note that we have not taken into account that the Roe variables permit a larger time step than the conservative variables.

Order of M W P = 2 P = 4 P = 8 P = 16

Time Roe var. general 29 32 44 107

Time cons. var. general 267 280 388 560

Time Roe var. const. e.vec. 6 7 12 29

Time cons. var. const. e.vec. 6 8 13 29

Table 1: Relative simulation time using conservative variables and Roe variables, respec-tively. One time unit is defined as the simulation time of a single deterministic problem with the same time-step as for the MW cases. Results are included for the codes designed for constant eigenvectors and for the same problem with the more general code which does not rely on the assumption of constant eigenvectors.

6. Conclusions

An intrusive formulation of the stochastic Euler equations based on Roe

variables is presented. A Roe average matrix for the standard

(31)

conditions stated by Roe under certain conditions.

The Legendre polynomial basis exactly represents the input uncertainty in our first test case, but it leads to oscillations around the discontinuity in stochastic space. The Haar wavelets do not represent the input uncertainty exactly in either test case, but are more robust to discontinuities.

The Roe variable formulation is robust for supersonic problems where the conservative variable formulation fails, but only for localized basis functions of the generalized chaos representation. For global Legendre polynomials, the discontinuities in stochastic space lead to oscillations and unphysical behavior of the solution and numerical breakdown. Haar wavelets are more robust in this respect, and do not yield oscillations around discontinuities in stochastic space. The robustness properties can be significantly improved with a stochastic basis that admit an eigenvalue decomposition with constant eigenvectors of the inner triple product matrix that occur frequently in the evaluation of pseudo-spectral operations. When this is the case, the Roe variables and conservative variables are similar in performance using the same time step.

For the general case where we do not assume an eigenvalue decomposition with constant eigenvectors, the Roe variable formulation leads to speedup compared to the conservative variable formulation. The relative speedup decreases with the order of generalized chaos since the total computational cost for high-order expansions is no longer dominated by spectral inversion and square root calculations. Instead, the main cost lies in the formation of spectral product matrices. However, for low order multi-wavelet expansions, the speedup is significant. The difference in computational time is mainly due to the pseudo-spectral operations of the numerical flux functions, especially if these are ill-conditioned.

We demonstrate the need for robust flux functions by presenting cases

where the standard MUSCL-Roe flux fails to capture the solution. The

design of a robust numerical method is also highly dependent on the choice of the stochastic basis. The Haar wavelets are not only more robust than Legendre polynomials for representation of discontinuities in stochastic space, but also admit the proof of existence of a Roe matrix and more specifically the hyperbolicity of the stochastic Galerkin formulation. This implies that the truncated problem mimics the original problem - a desirable feature.

If the representation of the initial function has not converged, the solu-tion at future times can not be accurate. The test case with uncertain initial shock location (case 2 in Section (5.2)) illustrates the need to find a

(32)

repre-sentation of uncertainty with fast decay of the coefficients of the generalized chaos expansion. An alternative to more accurate representation of the input uncertainty is to combine the intrusive Roe variable formulation presented here with multi-element methods, for instance in the manner presented in [9] or using adaptive methods [17].

Acknowledgements

This work is supported by King Abdullah University of Science and Tech-nology (KAUST).

Appendix A Generation of multi-wavelets

Algorithm 1 Generation of multi-wavelets

The goal is to generate the mother-wavelets W

i (ξ)}

Np

i=0 of (4). Start with a set of functions {f1

k} Np k=0, defined by fk1(ξ) =    ξk, ξ ∈ [−1, 0], −ξk, ξ ∈ [0, 1], 0, otherwise.

STEP 1: Orthogonalize with respect to the monomials 1, ..., ξNpusing

Gram-Schmidt to obtain {f2 k} Np k=0. STEP 2: for i← 0 to Np− 1 do

Make sure hfii+1ξNp+ii 6= 0 (otherwise reorder).

for j = i + 1 to N0 do w = hf i+2 j ξNp+ii hfi+2 i ξNp+ii

fji+3← fji+2− wfii+2 end for

end for

STEP 3: Orthogonalize {fii+2}Np

i=0 using G-S.

for i← Np to 0 do

ψW

i (ξ)← Apply Gram-Schmidt to fii+2.

end for

outputW

i (ξ)}

Np

(33)

Appendix B Proof of constant eigenvectors of A

Proposition 1. The matrix A defined by (7) for Haar wavelets j}Pj=1 has

constant eigenvectors for all P = 2Nr, N

r ∈ N.

Sketch of proof. We will use induction on the order P of wavelet chaos to show that the matrix A has constant eigenvectors for all orders P . In order to do this, we will need certain features of the structure of A. We can express A2P in terms of the matrix AP. Two properties of the triple producthψiψjψki will be used to prove that A indeed has the matrix structure presented. Property 1: Let i∈ {1, . . . , P }, j = k ∈ {P + 1, . . . , 2P } and let j0 and j00 be the progenies of j. Then

hψiψ2ji = hψiψj20i = hψiψj200i.

Property 2: Consider the indices i ∈ {P + 1, . . . , 2P }, j = k ∈ {2P +

1, . . . , 4P}. Then hψiψ2ji =    P1/2 if j first progeny of i −P1/2 if j second progeny of i 0 otherwise .

As the induction hypothesis, we assume that given AP for some P = 2Nr,

Nr ∈ N, the next order of triple product matrix A2P can be written

A2P =



AP QPMP

MPQTP Λ



where QP is the matrix of constant eigenvectors of AP satisfyingkQPk22 = P , MP = diag(wP +1, . . . , wP) and Λ is diagonal and contains the eigenvalues of

AP. Then, we have that

 AP QPMP MPQTP Λ   QP ±P1/2I  =  QPΛ± P1/2QMP P MP ± P1/2Λ  =  QP ±P1/2I  (Λ±P1/2M ),

so the eigenvalues and eigenvectors of A2P are given by Λ ± P1/2MP and

[QP,±P1/2I]T, respectively. For the next order of expansion, 4P , we have

A4P =      AP QPMP MPQTP Λ   QP ⊗ [1, 1] P1/2I ⊗ [1, −1]  M2P M2P  QP ⊗ [1, 1] P1/2I⊗ [1, −1] T Λ⊗ I2+ P1/2MP ⊗  1 0 0 −1      (B.1)

(34)

To see that this is indeed the structure of A4P, note that any non-zero matrix

entry not already present in A2P, can be deduced using properties 1 and 2,

and scaling the rows/columns by multiplication by the diagonal matrix M2P.

The structure of A4P follows from the construction of the Haar wavelet basis,

but we do not give a proof here.

One can verify that A4P given by (B.1) has the eigenvectors and

eigen-values Q4P =    QP ⊗ [1, 1] P1/2IP ⊗ [1, −1]  ±(2P )1/2I 2P   , Λ4P = Λ⊗ I2+ P1/2MP ⊗  1 0 0 −1  ± (2P )1/2M 2P,

so the eigenvectors are constant (but the eigenvalues are variable in the coef-ficients (wi)j through MP and M2P). The base cases P = 1, P = 2, can easily

be verified, so by induction AP˜ has constant eigenvectors for all P = 2Nr,

Nr ∈ N.

Appendix C Eigenvalue decompositions of A

C.1 Piecewise constant multi-wavelets (Haar wavelets)

C.1.1 Nr= 2 Q = 1 2     1 1 1 1 1 1 −1 −1 √ 2 √2 0 0 0 0 √2 √2    , Λ = diag     u0+ u1+ √ 2u2 u0+ u1− √ 2u2 u0− u1+ √ 2u3 u0− u1− √ 2u3    

(35)

C.1.2 Nr= 3 Q =                       1 1 1 1 1 1 1 1 1 1 1 1 −1 −1 −1 −1 √ 2 √2 √2 √2 0 0 0 0 0 0 0 0 √2 √2 √2 √2 2 −2 0 0 0 0 0 0 0 0 2 −2 0 0 0 0 0 0 0 0 2 −2 0 0 0 0 0 0 0 0 2 −2                       Λ = diag                       u0+ u1+ √ 2u2+ 2u4 u0+ u1 + √ 2u2− 2u4 u0+ u1 − √ 2u2+ 2u5 u0+ u1− √ 2u2− 2u5 u0− u1+ √ 2u3+ 2u6 u0− u1+ √ 2u3− 2u6 u0− u1− √ 2u3+ 2u7 u0 − u1 − √ 2u3− 2u7                      

(36)

C.2 Piecewise linear multi-wavelets C.2.1 Nr= 1 Q =           1 2 1 2 1 2 1 2 − √ 3+1 4 √ 3−1 4 √ 3+1 4 − √ 3−1 4 −1 2 1 2 − 1 2 1 2 − √ 3−1 4 − √ 3+1 4 √ 3−1 4 − √ 3+1 4           , Λ = diag           u0− √ 3+1 2 u1− u2− √ 3−1 2 u3 u0+ √ 3−1 2 u1+ u2− √ 3+1 2 u3 u0+ √ 3+1 2 u1− u2+ √ 3−1 2 u3 u0− √ 3−1 2 u1+ u2− √ 3+1 2 u3           C.2.2 Nr= 2 Q =                          1 √ 8 1 √ 8 1 √ 8 1 √ 8 1 √ 8 1 √ 8 1 √ 8 1 √ 8 √ 14+3√3 8 − √ 14+3√3 8 − √ 14−3√3 8 √ 14−3√3 8 √ 3+1 8√2 − √ 3+1 8√2 − √ 3−1 8√2 √ 3−1 8√2 − √ 3+1 4√2 − √ 3+1 4√2 − √ 3−1 4√2 − √ 3−1 4√2 √ 3−1 4√2 √ 3−1 4√2 √ 3+1 4√2 √ 3+1 4√2 √ 3+1 8√2 − √ 3+1 8√2 √ 3−1 8√2 − √ 3−1 8√2 − √ 14−5√3 8 √ 14−5√3 8 √ 14+5√3 8 − √ 14+5√3 8 0 −1 2 1 2 0 0 1 2 − 1 2 0 0 − √ 3−1 4 √ 3+1 4 0 0 − √ 3+1 4 √ 3−1 4 0 −1 2 0 0 1 2 1 2 0 0 − 1 2 √ 3−1 4 0 0 − √ 3+1 4 √ 3+1 4 0 0 − √ 3−1 4                         

(37)

Λ = diag                          u0 + q 14+3√3 8 u1− √ 3+1 2 u2+ √ 3+1 4 u3− √ 2u6+ √ 3−1 √ 2 u7 u0− q 14+3√3 8 u1− √ 3+1 2 u2− √ 3+1 4 u3− √ 2u4− √ 3−1 √ 2 u5 u0 − q 14−3√3 8 u1− √ 3−1 2 u2+ √ 3−1 4 u3+ √ 2u4+ √ 3+1 √ 2 u5 u0+ q 14−3√3 8 u1− √ 3−1 2 u2− √ 3−1 4 u3+ √ 2u6− √ 3+1 √ 2 u7 u0+ √ 3+1 4 u1+ √ 3−1 2 u2− q 14−5√3 8 u3+ √ 2u6+ √ 3+1 √ 2 u7 u0 − √ 3+1 4 u1+ √ 3−1 2 u2+ q 14−5√3 8 u3+ √ 2u4− √ 3+1 √ 2 u5 u0 − √ 3−1 4 u1+ √ 3+1 2 u2+ q 14+5√3 8 u3− √ 2u4+ √ 3−1 √ 2 u5 u0+ √ 3−1 4 u1+ √ 3+1 2 u2− q 14+5√3 8 u3− √ 2u6− √ 3−1 √ 2 u7                          [1] N. Wiener, The Homogeneous Chaos, American Journal of Mathematics

60(4) (1938) 897–936.

[2] R. G. Ghanem, P.D. Spanos, Stochastic finite elements: a spectral ap-proach, Springer-Verlag, New York, 1991.

[3] D. Xiu, G. E. Karniadakis, The Wiener–Askey polynomial chaos for stochastic differential equations, SIAM J. Sci. Comput. 24(2) (2002) 619–644.

[4] O. P. Le Maˆıtre, H. N. Najm, P. P. P´ebay, R. G. Ghanem, O. M. Knio,

Multi-resolution-analysis scheme for uncertainty quantification in chem-ical systems, SIAM J. Sci. Comput. 29 (2007) 864–889.

[5] T. Chantrasmi, A. Doostan, G. Iaccarino, Pad´e-Legendre approximants

for uncertainty analysis with discontinuous response surfaces, J. Com-put. Phys. 228 (2009) 7159–7180.

[6] O. P. Le Maˆıtre, O. M. Knio, H. N. Najm, R. G. Ghanem, Uncer-tainty propagation using Wiener-Haar expansions, J. Comput. Phys. 197 (2004) 28–57.

[7] O. P. Le Maˆıtre, H. N. Najm, R. G. Ghanem, O. M. Knio, Multi-resolution analysis of Wiener-type uncertainty propagation schemes, J. Comput. Phys. 197 (2004) 502–531.

(38)

[8] J. Tryoen, O. P. Le Maˆıtre, M. Ndjinga, A. Ern, Roe solver with entropy corrector for uncertain hyperbolic systems, J. Comput. Appl. Math. 235 (2010) 491–506.

[9] J. Tryoen, O. P. Le Maˆıtre, M. Ndjinga, A. Ern, Intrusive Galerkin methods with upwinding for uncertain nonlinear hyperbolic systems, J. Comput. Phys. 229 (2010) 6485–6511.

[10] P. Pettersson, G. Iaccarino, J. Nordstr¨om, Numerical analysis of the

Burgers’ equation in the presence of uncertainty, J. Comput. Phys. 228 (2009) 8394–8412.

[11] P. Pettersson, J. Nordstr¨om, G. Iaccarino, Boundary procedures for the

time-dependent Burgers’ equation under uncertainty, Acta Mathematica Scientia 30 (2010) 539–550.

[12] O. P. Le Maˆıtre, O. M. Knio, Spectral methods for uncertainty quan-tification, first edition, Springer, 2010.

[13] R. Abgrall, A simple, flexible and generic deterministic approach to uncertainty quantifications in non linear problems: application to fluid flow problems, Rapport de recherche, INRIA Bordeaux (2008).

[14] R. Abgrall, P. M. Congedo, C. Corre, S. Galera, A simple semi-intrusive method for uncertainty quantification of shocked flows, comparison with a non-intrusive polynomial chaos method, in: ECCOMAS CFD, Lis-bonne, Portugal, 2010.

[15] G. Po¨ette, B. Despr´es, D. Lucor, Uncertainty quantification for systems

of conservation laws, J. Comput. Phys. 228 (2009) 2443–2467.

[16] X. Wan, G. E. Karniadakis, Long-term behavior of polynomial chaos in stochastic flow simulations, Computer methods in applied mechanics and engineering, 195 (2006) 5582–5596.

[17] J. Tryoen, O. P. Le Maˆıtre, A. Ern, Adaptive anisotropic spectral stochastic methods for uncertain scalar conservation laws, SIAM J. Sci. Comput. 34(5) (2012) A2459–A2481.

[18] B. van Leer, Towards the ultimate conservative difference scheme. V - A second-order sequel to Godunov’s method, J. Comput. Phys. 32 (1979) 101–136.

(39)

[19] P. L. Roe, Approximate Riemann solvers, parameter vectors, and differ-ence schemes, J. Comput. Phys. 43 (1981) 357–372.

[20] M. K. Deb, I. M. Babuska, J. T. Oden, Solution of stochastic partial differential equations using Galerkin finite element techniques, Comput. Methods Appl. Mech. Eng. 190(48) (2001) 6359–6372.

[21] C. L. Pettit, P. S. Beran, Spectral and multiresolution Wiener expan-sions of oscillatory stochastic processes, J. Sound Vib. 294 (2006) 752– 779.

[22] X. Wan, G. E. Karniadakis, An adaptive multi-element generalized poly-nomial chaos method for stochastic differential equations, J. Comput. Phys. 209 (2005) 617–642.

[23] J. A. Witteveen, A. Loeven, H. Bijl, An adaptive Stochastic Finite Ele-ments approach based on Newton-Cotes quadrature in simplex eleEle-ments, Computers & Fluids 38(6) (2009) 1270-1288.

[24] B. K. Alpert, A class of bases in L2 for the sparse representations of integral operators, SIAM J. Math. Anal., 24 (1993) 246–262.

[25] B. J. Debusschere, H. N. Najm, P. P. P´ebay, O. M. Knio, R. G. Ghanem,

O. P. Le Maˆıtre, Numerical challenges in the use of polynomial chaos representations for stochastic processes, SIAM J. Sci. Comput. 26 (2005) 698–719.

[26] R. J. LeVeque, Finite-volume methods for hyperbolic problems, Cam-bridge University Press, CamCam-bridge, 2002.

[27] P. Pettersson, Q. Abbas, G. Iaccarino, J. Nordstr¨om, Efficiency of shock

capturing schemes for Burgers’ equation with boundary uncertainty, Enumath 2009, The eight European Conference on Numerical Math-ematics and Advanced Applications, June 29-July 3 (2009).

[28] M. J. D. Powell, A Fortran Subroutine for Solving Systems of Nonlinear Algebraic Equations, in: Numerical Methods for Nonlinear Algebraic Equations, P. Rabinowitz, 1970.

[29] P. J. Roache, Verification of Codes and Calculations, AIAA Journal 36(5) (1988) 696–702.

(40)

[30] L. Shunn, F. E. Ham, P. Moin, Verification of variable-density flow solvers using manufactured solutions, J. Comput. Physics 231(9) (2012) 3801–3827.

References

Related documents

undervisningen (Brumark, 2010). Genom att ha klassråd som den enda form av delaktighet för elever skapas det inte demokratiska medborgare av alla elever då många elever inte passar

We found that the risk of developing ARC and asthma before 10 years of age correlated with high SCORAD points during infancy. The SCORAD system has previously shown adequate

In the previous chapters I have shown how the separate diegetic narratives each correspond to different depths within the subconscious of the implied author, and also how

corresponding lateral variation in the c lattice-spacing. With temporal control of the substrate azimuthal orientation, nanospirals can be formed by, e.g., sequentially

In the talk and in the present extended abstract, we will first give a very concise introduction to stochastic partial differential equations with a particular focus on the

Although several studies have indicated similar molecular properties as being desirable for producing stable amorphous compounds (56,60,62,63,65), applying strict cut-off values

To clarify, using this transformation we know the possible realizations (the jump positions) of the random variable W ˜ n and we are back to the earlier discussed problem of finding

46 Konkreta exempel skulle kunna vara främjandeinsatser för affärsänglar/affärsängelnätverk, skapa arenor där aktörer från utbuds- och efterfrågesidan kan mötas eller