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
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)
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].
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
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.
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
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
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
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 ,
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 .
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
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.
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].
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.)
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)
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.
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
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
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].
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 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
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
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.
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
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
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
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.
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).
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
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
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
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
output{ψW
i (ξ)}
Np
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)
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
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
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
Λ = 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.
[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.
[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.
[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.