A posteriori error analysis and adaptive non-intrusive numerical schemes for systems of random conservation laws

This article considers one-dimensional random systems of hyperbolic conservation laws. Existence and uniqueness of random entropy admissible solutions for initial value problems of conservation laws, which involve random initial data and random flux functions, are established. Based on these results an a posteriori error analysis for a numerical approximation of the random entropy solution is presented. For the stochastic discretization, a non-intrusive approach, namely the Stochastic Collocation method is used. The spatio-temporal discretization relies on the Runge–Kutta Discontinuous Galerkin method. The a posteriori estimator is derived using smooth reconstructions of the discrete solution. Combined with the relative entropy stability framework this yields computable error bounds for the entire space-stochastic discretization error. The estimator admits a splitting into a stochastic and a deterministic (space-time) part, allowing for a novel residual-based space-stochastic adaptive mesh refinement algorithm. The scaling properties of the residuals are investigated and the efficiency of the proposed adaptive algorithms is illustrated in various numerical examples.


Introduction
Quantifying the influence of random model parameters as well as uncertain initial or boundary conditions has become an important task in computational science and engineering.Uncertainty Quantification (UQ) addresses this issue and provides a variety of mathematical methods to examine the influence of uncertain input parameters on numerical solutions and derived quantities of interest.
In this article we study (spatially) one-dimensional systems of random hyperbolic conservation laws, where the uncertainty stems from random initial data or from random flux functions.Based on Kružkov's work [23], a firm theory for random scalar conservation laws in several space dimensions has been developed [25,26,30].Compared to scalar equations, little is known about existence and uniqueness of entropy solutions for systems of hyperbolic conservation laws, especially in multiple space dimensions.For one-dimensional systems with initial data with small total variation, Glimm provided a proof for the global existence of entropy admissible solutions [17].Later, Bressan and coauthors [5,6,7] proved that the entropy admissible solutions constructed by the Glimm scheme (or equivalently by wave-front tracking) are unique in the sense that they are the only entropy admissible solutions satisfying additional stability properties such as certain bounds on the growth of their total variation.Based on these deterministic results, we prove the existence and uniqueness of so-called random entropy admissible solutions for one-dimensional systems of random hyperbolic conservation laws (Theorem 3.8).
Having established existence and uniqueness we approximate the random entropy admissible solution numerically.We discretize the random space by the Stochastic Collocation (SC) method [1,2,36].The method is non-intrusive, which means that the underlying deterministic solver does not need to be modified.Moreover, it is easily parallelizable and it avoids the problem of losing hyperbolicity for nonlinear hyperbolic systems, a major drawback of many intrusive methods, most notably the Stochastic Galerkin method [28].As specific deterministic solver we use the Runge-Kutta Discontinuous Galerkin method [11].
For nonlinear random hyperbolic conservation laws, discontinuities, both in physical and stochastic space, may appear in finite time.It is therefore favorable to locally increase the spatial and stochastic mesh resolution around the discontinuities in physical and stochastic space.Adaptive algorithms for the Stochastic Collocation method and for the Simplex Stochastic Collocation method have been considered in [19,35].The refinement criteria for these methods are based on heuristic considerations and are not immediately linked to the true numerical error.Moreover, they do not consider refinement in physical space.An approach which combines both, physical and stochastic refinement, has been introduced in [8], where the authors consider random boundary value problems for second order partial differential equations and use adjoint methods to derive separable error bounds for linear quantities of interest.They then use the corresponding residuals for local mesh refinement.For random elliptic problems the analysis of adaptive mesh refinements based on a posteriori error estimates is more advanced than for random hyperbolic conservation laws, cf.[18,31] and references therein.
In this work, we present as the first new contribution an a posteriori error analysis which is based on the following approach [24]: We view the numerical solution (or, more precisely, a reconstruction thereof) as the exact solution of a perturbed version of the original problem.The perturbation is given by a computable residual which acts as a source term.By using the appropriate stability theory for the problem at hand we can bound the difference between the exact and the numerical solution in terms of the residual.The suitable stability theory for systems of hyperbolic conservation laws is the relative entropy stability framework of Dafermos and DiPerna, see [12,Theorem 5.2.1 ].
The specific reconstructions, which we use, are based on reconstructions for deterministic problems suggested in [13], see also [16] for a modified reconstruction in terms of Stochastic Galerkin schemes.Using these reconstructions, we obtain a residual admitting a decomposition into a spatial and a stochastic part, which enables us to control the errors arising from spatial and stochastic discretization.Based on the a posteriori error estimate we exploit the residuals' structure and propose as the second novel contribution of this paper a residual-based space-stochastic adaptive numerical scheme.While for smooth solutions the estimator provides a reliable a posteriori error control, for discontinuous solutions it blows up under mesh refinement.However, the residuals precisely capture the positions of rarefaction waves, contact discontinuities and shocks.Therefore, our residual-based space-stochastic adaptive numerical scheme leads to a significant error reduction compared to uniform mesh refinement.Due to its non-intrusive structure our proposed method admits a straightforward parallelization, the residuals are on-the-fly computable and the resulting adaptive schemes efficiently decrease the numerical error compared to uniform mesh refinement.
This article is structured as follows: In Section 2 we describe our equation of interest.In Section 3 we first review the deterministic well-posedness theory for one-dimensional systems of hyperbolic conservation laws.We then introduce the notion of random entropy admissible solutions and establish existence and uniqueness under suitable assumptions on the random initial data and random flux function.Section 4 describes the stochastic discretization and we show how to obtain the reconstruction from our numerical solution.
In Section 5 we establish the a posteriori estimate and derive the splitting of the error estimator into a spatio-temporal and a stochastic part.Furthermore, we describe our space-stochastic adaptive numerical schemes.Finally, in Section 6, we provide various numerical examples for the Euler equations, where on the one hand we examine the scaling behavior of the corresponding residuals and on the other hand we compare the error reduction of our adaptive numerical algorithms to that of uniform mesh refinements and show that our adaptive schemes are indeed more efficient.

A Primer on Probability Theory
Let pΩ, F, Pq be a probability space, where Ω is the set of all elementary events ω P Ω, F is a σ-algebra on Ω and P is a probability measure.In the following we parametrize the uncertainty with a random vector ξ : Ω Ñ Ξ Ă R N with independent, absolutely continuous components, i.e. ξpωq " `ξ1 pωq, . . ., ξ N pωq ˘: Ω Ñ Ξ Ă R N .This means that for every random variable ξ i there exists a density function w i : R Ñ r0, 8q, such that ş R w i pyq dy " 1 and Prξ i P As " ş A w i pyq dy, for any A P BpRq, for all i " 1, . . ., N .Here BpRq is the Borel σ-algebra on R.Moreover, the joint density function w of the random vector ξ " pξ 1 , . . ., ξ N q can be written as wpyq " w i py i q @ y " py 1 , . . ., y N q J P Ξ.
The random vector induces a probability measure PpBq :" Ppξ ´1pBqq for all B P BpΞq on the measurable space pΞ, BpΞqq.This measure is called the law of ξ and in the following we work on the image probability space pΞ, BpΞq, Pq.
For a Banach space E and its Borel σ-algebra BpEq, we consider the weighted L p wspaces equipped with the norms }f } L p w pΞ;Eq :" ess sup yPΞ }f pyq} E , p " 8.

Statement of the Problem
We start with the following one-dimensional hyperbolic system of m P N (deterministic) conservation laws.
Here, upt, xq P U Ă R m is the vector of conserved quantities, F P C 2 pU; R m q is the flux function, and U Ă R m is the state space, which is assumed to be an open set and T P p0, 8q.We make the following assumptions on the initial condition and flux function.
(D1) The initial condition satisfies u 0 P L 1 pR; Uq.For our probabilistic equation of interest we admit in (IVP) random variations in the flux and initial datum, leading to the following one-dimensional system of m P N random conservation laws.
For the sake of simplicity we keep the same notation for the solution and for the flux as in (IVP) and make the following assumptions on the random initial condition and on the random flux function.Note that these assumption are the probabilistic versions of assumptions (D1) and (D2).
(R1) The uncertain initial condition satisfies u 0 P L 1 w pΞ; L 1 pR; Uqq.(R2) For almost every realization y P Ξ, the Jacobian D F p¨, yq has m distinct real eigenvalues, and each characteristic field is either linearly degenerate or genuinely nonlinear.Moreover, we assume that F P L 2 w pΞ; C 2 pU; R m qq.
3 Well-Posedness: Deterministic vs. Random Hyperbolic Conservation Laws In this section we first review some classical results for deterministic one-dimensional hyperbolic conservation laws.We then introduce the notion of a random entropy solution for (RIVP) and establish its existence and uniqueness based on the results for the deterministic hyperbolic conservation law (IVP).

Deterministic Hyperbolic Conservation Laws
Let us consider the deterministic initial value problem (IVP).We say that a strictly convex function η P C 2 pU; Rq and a function q P C 2 pU; Rq form an entropy/entropy-flux pair, if they satisfy D η D F " D q.We assume that the system (IVP) is endowed with at least one entropy/entropy-flux pair.We then define entropy admissible solutions in the following way.
Definition 3.1 (Entropy admissible solution).A function u P L 1 pp0, T q ˆR; Uq is called an entropy admissible solution of (IVP), if it satisfies the following conditions: 1.It is a weak solution, i.e.
for all φ P C 8 c pr0, T q ˆR; R m q.
2. It satisfies the weak entropy inequality: for all Φ P C 8 c pr0, T q ˆR; R `q.Theorem 3.2 ([6], Theorem 2).Provided (D2) holds, there exists a non-empty closed domain D Ă L 1 pR; Uq of integrable functions with small total variation and a semi-group Sptq : r0, 8q ˆD Ñ D, called Standard Riemann Semigroup (SRS), that is unique (up to its domain) and which has in particular the following properties: (i) There exists a constant L ą 0, such that for all u, v P D and for all s, t ě 0.
(ii) For u P D the function upt, xq :" pSptquqpxq is an entropy admissible solution of (IVP).It is the unique entropy admissible solution that is obtained as L 1 -limit of the wave-front tracking algorithm.Remark 3.3 (Uniqueness).While (IVP) may have several entropy admissible solutions there is one and only one entropy admissible solution induced by the SRS; in this sense entropy admissible solutions induced by SRS are unique.It was proven in [6] that the SRS-induced entropy admissible solution is the only entropy admissible solution satisfying certain additional stability properties, cf.[6, (A2),(A3)].Remark 3.4 (Domain of SRS).The domain of the SRS is discussed in [6, equation (1.3)].Note that it can always be replaced by a smaller set in order to make sure that additional conditions (such as (3.3)) hold.
Additionally, we will use the following result on the stability of the SRS.In particular, we can quantify how much the SRS-induced entropy admissible solution varies if the flux is changed.
for some suitable positive M P R and some compact set C Ă R m .For t ą 0 we denote by Spt, F q the SRS from Theorem 3.2, associated with the flux function F .Then there exists a constant C ą 0, such that for any flux function F , satisfying (D2) and Dp F q Ď DpF q, it holds that for all u P Dp F q.

Existence and Uniqueness of Random Entropy Solutions
We now consider the random hyperbolic conservation law (RIVP) and introduce the notion of a random entropy admissible solution for (RIVP).We say that η P L 1 w pΞ; C 2 pU; R m qq, q P L 1 w pΞ; C 2 pU; R m qq form a random entropy/entropy-flux pair if ηp¨, yq is strictly convex P-a.s.y P Ξ and if η and q satisfy D ηp¨, yq D F p¨, yq " D qp¨, yq, P-a.s.y P Ξ.We assume that the random conservation law (RIVP) is equipped with at least one random entropy/entropy-flux pair.
We define random entropy admissible solutions as path-wise (w.r.t.y) entropy admissible solutions of (RIVP).In this sense the notion of random entropy admissible solutions generalizes the notion of entropy admissible solutions in a similar way as the notion of random entropy solutions, introduced by Schwab and Mishra [26] generalizes the notion of entropy solutions according to Kružkov [23].

Definition 3.6 (Random entropy admissible solution). A function u P L 1
w pΞ; L 1 pp0, T qR ; Uqq is called a random entropy admissible solution of (RIVP), if it satisfies the following conditions: 1.It is a weak solution: P-a.s.y P Ξ and for all φ P C 8 c pr0, T q ˆR; R m q.
2. It satisfies the weak entropy inequality: P-a.s.y P Ξ and for all Φ P C 8 c pr0, T q ˆR; R `q.Remark 3.7.(i) Let u be a random entropy admissible solution, then for almost any fixed realization ỹ P Ξ, the function up¨, ¨, ỹq is an entropy admissible solution of the deterministic version of (RIVP) in the sense of Section 3.1.
(ii) The function upt, x, yq :" `Spt, F p¨, yqqu 0 p¨, yq ˘pxq, where tSpt, F p¨, yqqu tě0 is the SRS from Theorem 3.2 associated with the flux-function F p¨, yq, obviously satisfies (3.5) and (3.6) P-a.s.y P Ξ. Theorem 3.8 below discusses the regularity of u w.r.t.y, i.e. under which conditions u is indeed a random entropy admissible solution of (RIVP).
To ensure the existence of a random entropy admissible solution of (RIVP) by applying Theorem 3.2 and Theorem 3.5 path-wise in Ξ we make the following assumptions: (R3) We define D :" Ş yPΞ DpF p¨, yqq, where DpF p¨, yqq is the domain of the SRS from (3.3) in Theorem 3.5.We assume that D ‰ Ø and u 0 p¨, yq P D, P-a.s.y P Ξ.
Then u is a random entropy admissible solution of (RIVP).It is unique in the sense that it is the only random entropy admissible solution which path-wise coincides with the SRS-induced entropy admissible solution of the deterministic version of (RIVP).
Proof.The function u is path-wise the unique SRS-induced entropy solution of (RIVP) by construction.Note that we have assumed u 0 p¨, yq P D Ă Dp¨, F p¨, yqq, P-a.s.y P Ξ.It remains to show, that u is a random variable, i.e.
Hence, the mapping Sptq : pu, F q Q E 1 Ñ L 1 pR; R m q, Sptqpu, F q :" Spt, F qup¨q is continuous for all t ą 0. Due to the finite time horizon we immediately deduce that the mapping S : E 1 Ñ L 1 pp0, T q ˆR; R m q, Spu, F q :" Sp¨, F qup¨q is also continuous.Finally, it follows from assumptions (R1) and (R2) that the mapping S 0 : ´Ξ, BpΞq ¯Ñ ´E1 , BpE 1 q ¯, S 0 pyq :" pu 0 p¨, yq, F p¨, yqq is measurable.Thus, up¨, ¨, yq " Sp¨, F p¨, yqqu 0 p¨, yq " S ˝S0 pyq is a composition of measurable mappings and hence measurable itself.

Space-Time Stochastic Discretization and Reconstructions
A major goal of this paper is to prove an a posteriori estimate for a large class of numerical approximations of (RIVP).In particular, we consider numerical schemes that combine Stochastic Collocation (SC) with Runge-Kutta Discontinuous Galerkin (RKDG) schemes.To this end, we recapitulate the SC method for the (non-intrusive) discretization of the random space Ξ as for example in [1,36].Additionally, we recall the Multi-Element method which decomposes the random space Ξ into smaller elements to allow for a local interpolation in the random space [34].Finally, we describe a reconstruction of the numerical solution as a Lipschitz continuous function.The reconstruction will be used in the a posteriori error estimate in Section 5.

The Stochastic Collocation Method
The idea of the SC method is to approximate the random entropy admissible solution of (RIVP) by a polynomial interpolant in the random space, where the interpolant is supposed to satisfy (RIVP) at collocation points ty i u M i"0 Ă R, M P N. The exact solution up¨, ¨, y i q at a given collocation point y i , i " 0, . . ., M , is then approximated by a discrete solution u h p¨, ¨, y i q using any suitable numerical method.
For a multi-dimensional random space Ξ Ă R N , we define Ξ i :" ξ i pΩq and follow [2] to first define the space P q pΞq of tensor product polynomials of maximal degree q P N 0 by P q pΞq :" N â i"1 P q pΞ i q, P q pΞ i q :" tp : Ξ i Ñ R | p is a polynomial of degree qu.Remark 4.1.Our analysis does not depend on the structure of the approximating space, i.e. instead of considering a fixed polynomial degree q P N 0 for every random dimension we could also consider variable polynomial degrees q i P N, i " 1, . . ., N .Moreover, instead of using a tensor product space we could also consider the complete polynomial space, cf.[36] or sparse grids, cf.[9,27].
We let K :" tk " pk 1 , . . ., k N q J P N N 0 : k i ď q, i " 1 . . ., N u be the corresponding multi-index set and define the collocation points y k " py 1,k 1 , . . ., y N,k N q P Ξ, k P K.As a basis of P q pΞ i q we choose the Lagrange basis tl i,k u q k"0 associated with the collocation points ty i,k u q k"0 , such that l i,k py i,j q " δ k,j , @ j, k " 0, . . ., q, for all i " 1, . . ., N .We then define the multivariate Lagrange polynomials via l k py j q :" l 1,k 1 py 1,j 1 q ¨. . .¨lN,k N py N,j N q, j, k P K.
Using the collocation points ty k u kPK as input parameters in (RIVP) yields cardpKq " pq `1q N (uncoupled) collocated initial value problems: # B t upt, x, y k q `Bx F pupt, x, y k q, y k q " 0, pt, xq P p0, T q ˆΛ, up0, x, y k q " u 0 px, y k q, Λ P R . (CIVP) Here and in the following we consider Λ P tr0, 1s per , Ru.
Remark 4.2.The well-posedness result in Theorem 3.8 only covers Λ " R.However, the deterministic well-posedness results are based on local estimates and we, therefore, believe that it can be extended to cover the case Λ " r0, 1s per , as well.
Each of the deterministic hyperbolic systems in (CIVP) can be solved using the RKDG method described in Appendix A. For every collocation point y k we denote the corresponding numerical approximation at time t n " t n py k q by u n h p¨, y k q :" u h pt n , ¨, y k q.Let us assume that the time-partition tt n u Nt n"0 is the same for every collocation point ty k u kPK .The numerical approximation of (RIVP) at time t " t n can then be written as An important aspect of the SC method is the choice the collocation points ty i,k u q k"0 Ă Ξ i .Depending on the distribution of the random variable ξ i we choose the collocation points as zeros of the corresponding (orthogonal) chaos polynomials [37].For example, if ξ i " Upa, bq is uniformly distributed, we choose ty i,k u q k"0 to be the roots of the pq `1qth Legendre polynomial.For a Gaussian distribution we use the roots of the Hermite polynomials accordingly.
One can then approximate the mean of u n h px, ¨q via numerical quadrature, i.e.
Here w k are products of the corresponding one-dimensional weights.

The Multi-Element Stochastic Collocation Method
A major drawback of any global approximation approach in Ξ for hyperbolic conservation laws is that, due to the Gibbs phenomenon, the interpolant may oscillate for discontinuous solutions, cf.[28,33].To overcome this issue, we employ the Multi-Element (ME) approach as presented in [34], i.e. we subdivide Ξ into disjoint elements and consider a local approximation of (RIVP) on every random element.
For the ease of presentation we assume that Ξ " r0, 1s N , and let 0 " d 1 ă d 2 ă . . .ă d N Ξ `1 " 1 be a decomposition of r0, 1s.We define D n :" rd n , d n`1 q, for n " 1, . . ., N Ξ ´1, and D N Ξ :" rd N Ξ , d N Ξ `1s.Introducing the set M :" tm " pm 1 , . . ., m N q J P N N 0 : m i ď N Ξ , i " 1 . . ., N u allows us to define for m P M, the Multi-Element D m :" Hence, we consider a new local random variable ξ m : ξ ´1pD m q Ñ D m on the local probability space pξ ´1pD m q, F Xξ ´1pD m q, Pp¨|ξ ´1pD m qqq.Using Bayes' rule we can compute the local probability density function of ξ m via w ξ m :" w ξ py m |ξ ´1pD m qq " where Ppξ ´1pD m qq ą 0 for m P M can be assumed w.l.o.g., due to the independence of the corresponding random variables.We parametrize the uncertain input in (RIVP) using the local random variable ξ m and consider a local approximation on every D m at time t " t n , n " 0, . . ., N t , for all m P M. Here, ty m k u kPK Ă D m are the local collocation points in D m , m P M. The global Multi-Element Stochastic Collocation (ME-SC) approximation at time t " t n , can then be written as where χ Dm is the indicator function of D m .

Space-Time-Stochastic Reconstructions
For the space-time discretization of (CIVP) we use the RKDG framework from [10].For ease of presentation, we move the description of the RKDG scheme and the computation of its space-time reconstruction to Appendix A. As discussed in Appendix A, we have for each collocation point y k , k P K, a computable space-time reconstruction ûst py k q " ûst p¨, ¨, y k q P W 1 8 pp0, T q; V s p`1 X C 0 pΛqq of the numerical approximation u h p¨, ¨, y k q, where V s p`1 denotes the space of piece-wise polynomials of degree p `1 on a triangulation of Λ.This allows us to define the space-time residual as follows.
Definition 4.3 (Space-time residual).We call the function R st py k q :" R st p¨, ¨, y k q P L 2 pp0, T q ˆΛ; R m q, defined by R st pt, x, y k q :" B t ûst pt, x, y k q `Bx F pû st pt, x, y k q, y k q, (4.5) the space-time residual associated with the collocation point y k , for all k P K.
This residual is required in the subsequent analysis.In the next step we expand the space-time reconstruction into the corresponding random basis, i.e. in the Lagrange basis, to obtain the so-called space-time-stochastic reconstruction.Definition 4.4 (Space-time-stochastic reconstruction).We call the function ûsts P P q pΞqb `W 1 8 p0, T q; V s p`1 X C 0 pΛq ˘defined by ûsts pt, x, yq :" ÿ kPK ûst pt, x, y k ql k pyq, the space-time-stochastic reconstruction of the numerical approximation u h of (RIVP) (see (4.1)).
Similar to the space-time reconstruction, we may plug ûsts into the random conservation law (RIVP) to obtain the so called space-time-stochastic residual.Definition 4.5 (Space-time-stochastic residual).We define the space-time-stochastic residual R sts P L 2 w pΞ; L 2 pp0, T q ˆΛ; R m qq by R sts pt, x, yq :" B t ûsts pt, x, yq `Bx F pû sts pt, x, yq, yq.
We also need this residual for the upcoming error analysis.

A Posteriori Error Estimate and Adaptive Algorithms
As already mentioned in the introduction, our a posteriori error estimate relies on the relative entropy stability framework of Dafermos and DiPerna, see [12] and references therein.The relative entropy method allows to measure the L 2 -distance of two functions, one of them required to be Lipschitz continuous in space and time.This is the reason why we reconstructed the numerical solution and computed the space-time-stochastic reconstruction as a Lipschitz function.
Before stating the main a posteriori error estimate, we establish bounds on the derivatives of the flux function and the entropy, as they enter the upper bounds in the main estimate.Due to Assumption (R4) from Section 3.2 and the compactness of C, there exist P-a.s.y P Ξ constants 0 ă C F pyq ă 8 and 0 ă C η pyq ă C η pyq ă 8, such that, Here, for a generic function f , H u f denotes its Hessian matrix which contains all second order derivatives with respect to u.

A Posteriori Error Estimate and Error Splitting
We now have all ingredients together to state the following main a posteriori error estimate that can be directly derived from Theorem 5.5 in [15].Theorem 5.1 (A posteriori error bound for the numerical solution).Let u be a random entropy admissible solution of (RIVP).Let the reconstruction ûsts only take values in C.Then, the difference between u and the numerical solution u n h from (4.
In Theorem 5.1 the error between the numerical solution and the random entropy admissible solution is bounded by the error in the initial condition, the difference between the numerical solution and its reconstruction, and the contribution of the space-time stochastic residual R sts from Definition 4.5, quantified by E sts .However, we still cannot distinguish between errors that arise from discretizing the random space and the physical space.We, thus, are going to describe a splitting of the space-time-stochastic residual R sts into a spatial and a stochastic residual.Lemma 5.2 (Splitting of the space-time-stochastic residual).The space-time-stochastic residual R sts admits the decomposition with (5.4) R det and R stoch are called the deterministic and stochastic residual.
Proof.For every collocation point y k , k P K, we compute in Appendix A the space-time reconstruction ûst p¨, ¨, y k q which fulfills R st py k q " B t ûst py k q `Bx F pû st py k q, y k q. (5.
To simplify Theorem 5.1 let us assume that the eigenvalues of the Hessian H u ηpu, yq are bounded from above and below by positive numbers, for any u P C uniformly in Ξ.We let C η :" ess sup In our numerical examples, where we consider the compressible Euler equations, the dependence of the flux function F and η on y is explicitly known and we can compute the constants C η , C η , C F numerically.
The following corollary is a simple consequence of the splitting in Lemma 5.2.

Corollary 5.3 (A posteriori error bound with error splitting and simplified bounds).
Let u be a random entropy admissible solution of (RIVP).Then, the difference between u and the numerical solution u n h from (4.1) satisfies }upt n , ¨, ¨q ´un h p¨, ¨q} 2 for n " 0, . . ., N t and for all P-measurable sets Ξ Ď Ξ.Here, w p Ξ;L 2 pp0,tnqˆΛq , (5.9) E sts 0 :" }u 0 p¨, ¨q ´û sts p0, ¨, ¨q} 2 L 2 w p Ξ;Λq . (5.11) Remark 5.4.(i) The residual R det in (5.4) interpolates spatio-temporal residuals and contains information about the discretization error in physical space, i.e. the spacetime resolution of (CIVP) using the RKDG method.In contrast to R det , the stochastic residual R stoch in (5.5) indicates the quality of the interpolation in stochastic space.
(ii) In order for the upcoming space-stochastic adaptive algorithm based on E det , E stoch to be efficient, we need E det to depend solely on the space-time discretization and to be independent of the stochastic discretization.Similarly, we need E stoch to decay when the stochastic resolution is increased but to be independent of the space-time discretization.In Remark 5.5 we prove that E det is indeed unaffected by stochastic refinement.Our numerical experiment in Section 6.3.1 also shows that E stoch is unaffected by spatial refinement.
(iii) The scaling properties of E det , resp.R st py k q, were studied in [13].Currently we are not able to prove any of the scaling properties of E stoch w.r.t. to q and the number of Multi-Elements.However, our numerical experiments show that E stoch scales as desired, i.e.E stoch shows the same qualitative behavior as the stochastic interpolation error of the exact solution.
(iv) As described in [13] and [16], R det scales with 1 h in the vicinity of shocks and contact discontinuities, i.e., it blows up under spatial mesh refinement in these areas, although the numerical solution converges towards the exact solution.Hence, we only have reliable a posteriori error control for smooth solutions of (RIVP).However, as R det precisely captures the positions of rarefaction waves, contact discontinuities and shocks we use R det and R stoch , resp.E det and E stoch , as local indicators for our adaptive mesh refinement algorithms described in Section 5.2.
Remark 5.5 (Uniformity of the deterministic residual in Ξ).As noted above, the collocation points y k are chosen to be the zeros of the corresponding orthogonal polynomial depending on the distribution of ξ.The deterministic residual R det from (5.4) consists of Lagrange polynomials associated with the corresponding collocation points, thus Gaussian quadrature in Ξ yields Hence, E det inherits the convergence order of R st py k q.

Adaptive Algorithms
The splitting of the space-time-stochastic residual into a deterministic and a stochastic residual helps us in developing adaptive numerical schemes where we use the residuals as local error indicators for spatial and stochastic mesh refinement.We describe the deterministic spatially adaptive algorithm, which we use to solve (CIVP) for every collocation point y k , k P K. We slightly abuse the notation from (5.9) and write for every physical cell I P T n , E det k pt n , t n`1 , Iq :" }R st py k q} L 2 pptn,t n`1 qˆIq , which is the cell-wise indicator for the spatial refinement in Λ.
The local physical mesh refinement is achieved by uniformly dividing one cell into two new children cells or merging two cells into one parent cell.To mark elements for refinement we compute the deterministic residual E det k pt n , t n`1 , Iq on every cell I P T n and based on the residual we mark a fixed fraction of the cells for refinement.To coarsen the mesh, we can only merge cells that have the same parent element and both siblings are marked for coarsening.For coarsening we also choose a fixed fraction of all elements according to the local residual E det k pt n , t n`1 , Iq, cf.[20].Additionally, each cell is augmented with a variable denoting its current mesh-level which is initially zero.We fix a maximum mesh-level L P N, to restrict the fineness of the adaptive mesh.The algorithm reads as follows: Algorithm 1 Deterministic h-adaptive algorithm Input: final time T , max mesh-level L, initial mesh T 0 1: Compute u n`1 h on the current mesh T n using Algorithm 3 (see Appendix A) , Iq for I P T n and mark a fixed fraction of the elements for refinement and coarsening a: Refinement: If the cell's mesh-level is L do nothing.Else divide it uniformly into two new cells and increase the two new cells' mesh-level by one b: Coarsening: If the cell's mesh-level is zero do nothing.Else check if its sibling is marked for coarsening.If yes merge the two cells into one and decrease its mesh-level by one 3: Project u n`1 h onto the new mesh T n`1 using the L 2 -projection 4: If t n`1 ă T go to step 1 Remark 5.6.After every projection step in line three of Algorithm 1 we apply the TVBM slope limiter ΛΠ h from Appendix A.
In the numerical experiments in Section 6 we observed that setting the refinement and coarsening fractions to 1% and 20% provided the best error reduction.The finest mesh level will be L " 3. Next, we describe the stochastic adaptive algorithm for the ME-SC method from Section 4.2 using the stochastic residual E stoch as local indicator for stochastic refinement.
Algorithm 2 Stochastic N Ξ -Adaptive Algorithm Input: initial number of Multi-Elements M Ξ , max no. of Multi-Elements N Ξ , q `1 number of collocation points in each stochastic dimension 1: For every Multi-Element D m compute pq `1q N numerical samples using Algorithm 1 2: Compute E stoch pT q on every Multi-Element D m and uniformly subdivide the Multi-Element with the biggest residual, set M Ξ :" M Ξ `p2 N ´1q 3: If M Ξ ă N Ξ compute M samples on every new Multi-Element and go to 2

Numerical Examples
In this section we present various numerical examples concerning the scaling properties of the residuals and the performance of the adaptive algorithms.In Section 6.1 and Section 6.3 we examine the scaling properties of E det and E stoch .Sections 6.2, 6.4 and 6.5 assess the efficiency of our proposed adaptive algorithms.
As numerical solver we use the RKDG Code Flexi [21].The DG polynomial degrees will always be one or two and for the time-stepping we use the low storage SSP RKmethod of order three as in [22].The time-reconstruction is also of order three.As numerical fluxes we choose either the Lax-Wendroff numerical flux Gpu, vq :" F pwpu, vqq, wpu, vq :" or the Lax-Friedrichs numerical flux Gpu, vq :" 1 2 ´F puq `F pvq ¯`λpv ´uq.( In our example, the uncertainty is uniformly distributed.Therefore, we use the zeros of the Gauß-Legendre polynomials as collocation points.Computing E det , E stoch requires computing integrals, we approximate them via Gauß-Legendre quadrature where we use seven points in time, ten points in physical space and ten points in random space, except for Example 6.3, where for the global interpolation the number of quadrature points in random space will be 2q. In the following experiments we consider as instance of (RIVP) the one-dimensional compressible Euler equations for the flow of an ideal gas, which are given by where ρ describes the mass density, m the momentum and E the energy of the gas.The constitutive law for pressure p reads with the adiabatic constant γ " 1.4 if not specified otherwise.In the following figures we refer to the quantity }mpT, ¨, ¨q ´mNt h p¨, ¨q} L 2 w pΞ;L 2 pΛqq at final computational time T as numerical error, unless otherwise stated.We also plot the residuals E det pT q and E stoch pT q as in (5.9) and (5.10) from the momentum equation.Remark 6.1.Due to the structure of the flux Jacobian for the Euler equations (6.3), the first component of the stochastic residual R stoch from (5.5) vanishes when considering the Euler equations without source term.We therefore use the residuals for the momentum and the energy balance as indicators for our space-stochastic mesh refinements.

A Deterministic Problem with Smooth Solution
In this numerical example, we study the scaling properties of the deterministic residual E det from (5.9) for a uniform spatial mesh refinement.To this end, we construct a smooth exact solution by introducing an additional source term into the Euler equations.The exact solution reads as follows ¨ρpt, x, yq mpt, x, yq Ept, x, yq The numerical solution is computed on the domain Λ " r0, 1s per up to T " 0.5 and we use the Lax-Wendroff numerical flux (6.1).
In Table 1 we present the numerical error and the residual E det from (5.9) for the smooth solution (6.4) for DG polynomial degrees one and two.We can see that the error and the residual converge with the correct order of convergence, which is p `1, where p is the DG polynomial degree.
p   The numerical solution is computed on the domain Λ " r0, 1s up to T " 0.2 using the Lax-Friedrichs flux (6.2) and a DG polynomial degree of two.In this example we use exact boundary conditions.In Figure (1a) we compare the L 1 pΛq-and L 2 pΛq-error at time T between the numerical solution and the exact solution obtained with an exact Riemann solver [3].We can see that for the same number of spatial cells N s , the numerical error obtained with the adaptive numerical algorithm is smaller than for the uniform mesh refinement.The adaptive algorithm is also computationally more efficient than the uniform algorithm, which can be seen in the error vs. cpu time plot in Figure (1b).Remark 6.2.As discussed in Remark 5.4 (iv), R det scales with 1 h in the vicinity of shocks and contact discontinuities, i.e., it blows up under spatial mesh refinement in these areas.Thus, if we view the residual as an error indicator, it severely over-estimates the error so that it is to be called "inefficient" in these areas, according to the nomenclature of e.g.[32].From the point of view of mesh adaptation however, refinement based on R det leads to a reasonable refinement strategy that yields a considerable improvement in error decay compared to uniform mesh refinement (cf. Figure (1)).In particular, R det precisely captures the positions of rarefaction waves, contact discontinuities and shocks.
Over-estimating the error at discontinuities leads to maximal refinement at discontinuities and some refinement strategies for hyperbolic conservation laws suggest a maximal refinement close to shocks [29].

A Stochastic Problem with Smooth Solution
In this section we focus on the scaling properties of the stochastic residual for a one-and two-dimensional random space Ξ and a random flux function.

A One-Dimensional Random Space, q-Refinement
We modify the exact solution from Section 6.1 in the following way, ¨ρpt, x, yq mpt, x, yq Ept, x, yq The numerical solution is computed on Λ " r0, 1s per up to T " 0.2, the uncertainty y stems from an uniform distribution, i.e. ξ " Up0, 8q.We consider two different spatial meshes consisting of N s " 32 and 512 elements respectively, a DG polynomial degree of p " 2 and we use the Lax-Wendroff numerical flux (6.1).In this numerical example we globally approximate the function (6.5) in Ξ, i.e. we increase the polynomial degree q and consider one ME.Figure (2) shows the behavior of the error and the spatial, resp.stochastic residual, when we globally interpolate the smooth function (6.5).We see that the stochastic residual E stoch exhibits spectral convergence.Also the numerical error exhibits spectral convergence until it starts to stagnate because of the spatial resolution error.This is the correct behavior of the stochastic residual as we are globally increasing the polynomial degree in the random space and, therefore, expect spectral convergence with increasing  polynomial degree.We also observe that the exponential convergence of E stoch is not altered by a finer or coarser space discretization.Moreover, the deterministic residual E det is unaffected by the increasing resolution in the random space, which we expect from the residual's splitting into a space-time and a stochastic part.

Mesh Refinement in Ξ and Random Flux Function
In this example we examine the scaling properties of E stoch under mesh refinements for a two-dimensional random space Ξ Ă R 2 .We consider the same smooth function as in Section 6.3.1, ¨ρpt, x, y 1 q mpt, x, y 1 q Ept, x, y 1 q with ξ 1 " Up0, 8q.Moreover, we consider a random adiabatic constant.We assume that γ " ξ 2 " Up1.4,1.6q and thus the flux function is also random.The randomness of the adiabatic-constant corresponds to considering a gas mixture of uncertain composition.
The numerical solution is computed on Λ " r0, 1s per up to T " 0.2.We consider a fixed spatial mesh consisting of N s " 32 elements.For the ME-SC method we perform a linear and a quadratic interpolation, i.e. q P t1, 2u. Figure (3) illustrates the behavior of the stochastic residual E stoch , when we consider a local interpolation, i.e., when we consider the ME method from Section 4.2.We observe  that for a local linear and quadratic interpolation, i.e. q P t1, 2u, the stochastic residual converges approximately with the expected rate of convergence, which is pq`1q{2, cf.[34].
Like for the q-refinement in Section 6.3.1, the deterministic residual E det stays constant, when we increase the number of MEs.

Stochastic Adaptivity: Stochastic Problem with Discontinuous Solution
We apply the stochastic adaptive Algorithm 2 without spatial adaptivity to a solution which has a discontinuity in the random variable and compare the results with uniform space-stochastic mesh refinements.We therefore consider the following discontinuous function, ¨ρpt, x, y 1 , y 2 q mpt, x, y 1 , y 2 q Ept, x, y 1 , y 2 q ‚" ¨1 `Apy 1 , y 2 q cosp4πpx ´y1 tqq is a discontinuous amplitude.For the spatial domain Λ " r0, 1s per we use N s " 32 elements and a DG polynomial degree of two.The solution is computed up to T " 0.2 using the Lax-Wendroff numerical flux (6.1) and for the uncertainty we assume that ξ 1 , ξ 2 " Up0, 1q.For the ME-SC method we consider a linear interpolant, i.e. q " 1.

Space-Stochastic Adaptivity: An Uncertain Riemann Problem
Finally, we assess the efficiency of the space-stochastic adaptive algorithm by considering a random Riemann Problem.The initial data for this problem reads as follows ρpt " 0, x, yq " 1 mpt " 0, x, yq " # y 1 , x ď 0.5 y 2 , x ą 0.5 ppt " 0, x, yq " 1, where ξ 1 , ξ 2 " Up´1, 1q and Λ " r0, 1s.We compare the space-stochastic adaptive Algorithm 2 with uniform refinement, both in physical and random space.For this problem we use the Lax-Friedrichs numerical flux (6.2) and for the uniform spatial mesh we consider N s " 512 spatial elements.As for the Sod Shock Tube problem in Section 6.2 we prescribe exact boundary conditions.For the adaptive algorithm we always start on a spatial mesh consisting of N s " 256 elements.The DG polynomial degree is two and we consider a linear interpolation in the random space, i.e. q " 1.The solution is computed up to T " 0.2.The error is measured in the expected value rather than the L 2 w pΞ; L 2 pΛqq-norm.Note that we do not have an exact solution at hand for this problem, but due to Jensen's inequality, } EpupT, ¨, ¨qq ´Epu Nt h p¨, ¨qq} w pΞ;L 2 pΛqq .The reference expectation EpupT, ¨, ¨qq is computed using a Monte-Carlo method with an exact Riemann solver with 500000 samples.In Figure (5a) we show the numerical error as in (6.7) and we also consider the error } EpupT, ¨, ¨qq ´Epu Nt h p¨, ¨qq} L 1 pΛq for an increasing number of MEs, i.e. for increasing N Ξ .We can see that the adaptive algorithm decreases the error considerably faster than the uniform refinement.This is also depicted in the cpu time vs. error plot (Figure (5b)), where we can see that the adaptive algorithm reaches an absolute error in significantly less computational time than the uniform algorithm.This demonstrates, in particular, the efficiency of our proposed method.identify x 0 " x Ns to account for periodic boundary conditions.Further let 0 " t 0 ă t 1 ă . . .ă t Nt " T be a temporal decomposition of r0, T s and define ∆t n :" pt n`1 ´tn q, ∆t " max n ∆t n .With each time-interval pt n , t n`1 s we associate a (possibly different) partition T n and associated DG space V s p,n :" tv : Λ Ñ R m | v | I P P p pI, R m q, for all I P T n u.With L V s p,n we denote the L 2 -projection mapping into the DG space V s p,n .Following [14] we call the function u h a generalized semi-discrete DG approximation of (CIVP) if it satisfies for u ´1 h :" L V s p,0 u 0 the following equations.For every n " 0, . . ., N t , u n h | rtn,t n`1 s P C 1 ppt n , t n`1 q; V s p,n q X C 0 prt n , t n`1 s; V s p,n q, u n h pt n q " L V s p,n u n´1 h pt n q, Ns´1 ÿ i"0 x i`1 L n h pu n h q ¨ψh dx @ψ h P V s p,n , (DG) where L n h : V s p,n Ñ V s p,n is defined by Gpvpx í q, vpx ì qq ¨rr ψ h ss i , @v, ψ h P V s p,n .(A.1) The numerical solution u h is defined through u h p0q :" u ´1 h and u h | ptn,t n`1 s :" u n h | ptn,t n`1 s .Here, G : R m ˆRm Ñ R m denotes a numerical flux, the spatial traces are defined as ψpx ˘q :" lim hOE0 ψpx ˘hq and rr ψ h ss i :" pψ h px í q ´ψh px ì qq are jumps.
The initial-value problem (DG) can now be solved numerically by any single-or multistep method.We focus on K-th order Runge-Kutta time-step methods as in [11,22].Furthermore, ΛΠ h : R m Ñ R m is the TVBM minmod slope limiter from [11].Then, the complete S-stage time-marching algorithm for given n-th time-iterate u n h pt n q P V s p,n reads as follows: Algorithm 3 TVBM Runge-Kutta Time-Step 1: Set u p0q h = u n h pt n q. 2: for j " 1, . . ., S do `βjl α jl ∆t n L n h pu plq h q. 4: end for 5: Set u n h pt n`1 q " u pSq h .
We define the spatial reconstruction which is applied to the temporal reconstruction ût pt, ¨q for each t P p0, T q using the function w (cf.[13,15]).
Definition A.2 (Space-time reconstruction).Let ût be the temporal reconstruction of a sequence tu n h u Nt n"0 of solutions of the fully discrete scheme of (DG) using a numerical flux satisfying (A.2) or (A.3).The space-time reconstruction ûst pt, ¨q P V s p`1 is defined as the solution of Ns´1 ÿ i"0 x i`1 ż x i pû st pt, ¨q ´û t pt, ¨qq ¨ψ dx " 0 @ψ P V s p´1 , ûst pt, x k q " wpû t pt, x ḱ q, ût pt, x k qq @ k " 0, . . ., N s .
We have the following property of the space-time reconstruction.
cpu time: uniform vs. adaptive

Figure 1 :
Figure 1: Error plot for the deterministic Sod shock tube problem.Example 6.2

Figure 2 :
Figure 2: Error plot for stochastic smooth problem.Example 6.3.1

Figure 4 :
Figure 4: Error plot for discontinuous stochastic problem.Example 6.4 In Figure (4a) we plot the error and the spatial resp.stochastic residual versus the number of MEs and in Figure (4b) we show the error of the uniform and adaptive method versus cpu time.In Figure (4a) we can observe that for the uniform stochastic refinement, both the error and the stochastic residual E stoch converge with a rate of approximately 1{4.This is in accordance with what we expect when interpolating a two-dimensional discontinuous function.For the adaptive refinement the error and the residual exhibit a rate of convergence of approximately 1{2.The advantage of the stochastic adaptive algorithm is also reflected in Figure (4b), where we reach an error reduction in significantly less time compared to uniform refinement.

α
jl w jl h ¯, w jl h " u plq h