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

TL;DR

This paper develops an a posteriori error analysis and adaptive non-intrusive numerical schemes for one-dimensional random hyperbolic conservation laws, enabling efficient and reliable simulations with stochastic and spatial adaptivity.

Contribution

It introduces a novel residual-based adaptive mesh refinement algorithm for stochastic conservation laws using the Stochastic Collocation method and Runge--Kutta Discontinuous Galerkin discretization.

Findings

The error estimator effectively separates stochastic and deterministic errors.

The adaptive algorithm improves computational efficiency in numerical experiments.

Numerical results demonstrate the estimator's accuracy and the method's robustness.

Abstract

In this article we consider one-dimensional random systems of hyperbolic conservation laws. We first establish existence and uniqueness of random entropy admissible solutions for initial value problems of conservation laws which involve random initial data and random flux functions. Based on these results we present an a posteriori error analysis for a numerical approximation of the random entropy admissible solution. For the stochastic discretization, we consider a non-intrusive approach, the Stochastic Collocation method. The spatio-temporal discretization relies on the Runge--Kutta Discontinuous Galerkin method. We derive the a posteriori estimator using continuous 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. We conclude with various numerical examples investigating the scaling properties of the residuals and illustrating the efficiency of the proposed adaptive algorithm.

Tables

Table 1. Table 1: L 2 superscript 𝐿 2 L^{2} -error, residual and experimental order of convergence (eoc). Example 6.1 .

$p=1$
$N_s$	error	eoc	${ℰ}^{det}\left(T\right)$	eoc
16	1.2771e-02	-	1.1821e-01	-
32	4.5795e-03	1.48	4.5907e-02	1.36
64	1.4193e-03	1.69	1.5947e-02	1.53
128	3.5089e-04	2.02	4.3912e-03	1.86
256	7.9658e-05	2.14	1.3655e-03	1.69
512	1.9625e-05	2.02	3.7220e-04	1.88
1024	4.8984e-06	2.00	1.0804e-04	1.78

Equations

$$w(y)=\prod _{i=1}^N{w}_i(y_i)\forall y=(y_1,\dots ,y_N{)}^{\top }\in \Xi .$$

$$\displaystyle\|f\|_{L_{w}^{p}(\Xi;E)}:=\begin{cases}\Big{(}\int\limits_{\Xi}\|f(y)\|_{E}^{p}~{}w(y)\mathrm{d}y\Big{)}^{1/p}=\operatorname{\mathbb{E}}\Big{(}\|f\|_{E}^{p}\Big{)}^{1/p},\quad&1\leq p<\infty\\ \operatorname{ess~{}sup}\limits_{y\in\Xi}\|f(y)\|_{E},&p=\infty.\end{cases}$$

$$\{\begin{array}{ll}{\partial }_{t}u(t,x)+{\partial }_{x}F(u(t,x))=0,& (t,x)\in (0,T)\times R,\\ u(0,x)=u^0(x),& x\in R.\end{array}$$

$$\{\begin{array}{ll}{\partial }_{t}u(t,x,y)+{\partial }_{x}F(u(t,x,y),y)=0,& (t,x,y)\in (0,T)\times R\times \Xi ,\\ u(0,x,y)=u^0(x,y),& (x,y)\in R\times \Xi .\end{array}$$

$$\displaystyle\int\limits_{0}^{T}\int\limits_{\operatorname{\mathbb{R}}}\Big{(}u(t,x)\cdot\partial_{t}\phi(t,x)$$

$$=-\underset{R}{\int }{u}^0(x)\cdot \varphi (0,x)\mathrm{d}x,$$

$$\displaystyle\int\limits_{0}^{T}\int\limits_{\operatorname{\mathbb{R}}}\Big{(}\eta(u(t,x))\partial_{t}\Phi(t,x)$$

$$\ge -\underset{R}{\int }\eta (u^0(x))\Phi (0,x)\mathrm{d}x,$$

$$\Vert S(s)\overline{u}-$$

$$D(F)\subseteq \{u\in L^1(R,C)|TV(u)\le M\},$$

$$\Vert S(t,F)u-S(t,\tilde{F})u{\Vert }_{L^1(R,R^m)}\le Ct\Vert DF-D\tilde{F}{\Vert }_{C^0(U,R^m)},$$

$$\displaystyle\int\limits_{0}^{T}\int\limits_{\operatorname{\mathbb{R}}}\Big{(}u(t,x,y)\cdot\partial_{t}\phi(t,x)$$

$$=-\underset{R}{\int }{u}^0(x,y)\cdot \varphi (0,x)\mathrm{d}x,$$

$$\displaystyle\int\limits_{0}^{T}\int\limits_{\operatorname{\mathbb{R}}}\Big{(}\eta(u(t,x,y),y)\partial_{t}\Phi(t,x)$$

$$\ge -\underset{R}{\int }\eta (u^0(x,y),y)\Phi (0,x)\mathrm{d}x,$$

$$\Big{(}\Xi,\mathcal{B}(\Xi)\Big{)}\ni y\mapsto u(\cdot,\cdot,y)\in\Big{(}L^{1}((0,T)\times\operatorname{\mathbb{R}};\operatorname{\mathbb{R}}^{m}),\mathcal{B}((L^{1}((0,T)\times\operatorname{\mathbb{R}};\operatorname{\mathbb{R}}^{m})))\Big{)}$$

$$E_1:=L^1(R;R^m)\times C^2(U;R^m),$$

$$\Vert (\overline{u},\overline{F}){\Vert }_{E_1}:=\Vert \overline{u}{\Vert }_{L^1(R;R^m)}+\Vert \overline{F}{\Vert }_{C^2(U;R^m)}.$$

$$\Vert S$$

$$\le$$

$$+\Vert S(t,F(\cdot ,\tilde{y}))u^0(\cdot ,y)-S(t,F(\cdot ,\tilde{y}))u^0(\cdot ,\tilde{y}){\Vert }_{L^1(R;R^m)}$$

$$\le$$

$$\mathrm{S}:E_1\to L^1((0,T)\times R;R^m),\mathrm{S}(\overline{u},\overline{}F):=S(\cdot ,\overline{F})\overline{u}(\cdot )$$

$$\mathrm{S}_{0}:\Big{(}\Xi,\mathcal{B}(\Xi)\Big{)}\to\Big{(}E_{1},\mathcal{B}(E_{1})\Big{)},~{}\mathrm{S}_{0}(y):=(u^{0}(\cdot,y),F(\cdot,y))$$

$${\mathcal{P}}_q(\Xi ):=\underset{i=1}{\overset{N}{⨂}}{\mathcal{P}}_q({\Xi }_i),{\mathcal{P}}_q({\Xi }_i):=\{p:{\Xi }_i\to R|p\text{ is a polynomial of degree }q\}.$$

$$l_{i,k}(y_{i,j})={\delta }_{k,j},\forall j,k=0,\dots ,q,$$

$$l_{\mathrm{k}}(y_{\mathrm{j}}):=l_{1,k_1}(y_{1,j_1})\cdot \dots \cdot l_{N,k_N}(y_{N,j_N}),\mathrm{j},\mathrm{k}\in \mathcal{K}.$$

$$\{\begin{array}{ll}{\partial }_{t}u(t,x,y_{\mathrm{k}})+{\partial }_{x}F(u(t,x,y_{\mathrm{k}}),y_{\mathrm{k}})=0,& (t,x)\in (0,T)\times \Lambda ,\\ u(0,x,y_{\mathrm{k}})=u^0(x,y_{\mathrm{k}}),& \Lambda \in R.\end{array}$$

$$u_h^n(x,y):=\sum _{\mathrm{k}\in K}{u}_h^n(x,y_{\mathrm{k}})l_{\mathrm{k}}(y).$$

$$\displaystyle\operatorname{\mathbb{E}}\Big{(}u_{h}^{n}(x,\cdot)\Big{)}=\int\limits_{\Xi}u_{h}^{n}(x,y)~{}p_{\xi}(y)\mathrm{d}y\approx\sum\limits_{\mathrm{k}\in\operatorname{\mathcal{K}}}u_{h}^{n}(x,y_{\mathrm{k}})w_{\mathrm{k}}.$$

$$w_{{\xi }^{\mathrm{m}}}:=w_{\xi }(y^{\mathrm{m}}|{\xi }^{-1}(D_{\mathrm{m}}))=\frac{w_{\xi }(y^{\mathrm{m}})}{\mathbb{P}({\xi }^{-1}(D_{\mathrm{m}}))},y^{\mathrm{m}}\in D_{\mathrm{m}},$$

$$u_h^{\mathrm{m}}(t_n,x,y)=\sum _{\mathrm{k}\in K}{u}_h(t_n,x,y_{\mathrm{k}}^{\mathrm{m}})l_{\mathrm{k}}^{\mathrm{m}}(y),$$

$$u_h(t_n,x,y)=\sum _{\mathrm{m}\in \mathcal{M}}{u}_h^{\mathrm{m}}(t_n,x,y){\chi }_{D_{\mathrm{m}}}(y)=\sum _{\mathrm{m}\in \mathcal{M}}\sum _{\mathrm{k}\in K}{u}_h(t_n,x,y_{\mathrm{k}}^{\mathrm{m}})l_{\mathrm{k}}^{\mathrm{m}}(y){\chi }_{D_{\mathrm{m}}}(y),$$

Full text

Key words: hyperbolic conservation laws, uncertainty quantification, a posteriori error estimates, stochastic collocation method, discontinuous Galerkin method, adaptive mesh refinement

AMS subject classifications: Primary 35L65, 35R60; Secondary 65M15, 65M60, 65M70

1 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, 30, 26]. 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 [6, 7, 5] 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, 36, 2]. 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.

2 Preliminaries and Notation

2.1 A Primer on Probability Theory

Let $(\Omega,\operatorname{\mathcal{F}},\mathbb{P})$ be a probability space, where $\Omega$ is the set of all elementary events $\omega\in\Omega$, $\operatorname{\mathcal{F}}$ is a $\sigma$-algebra on $\Omega$ and $\mathbb{P}$ is a probability measure. In the following we parametrize the uncertainty with a random vector $\xi:\Omega\to\Xi\subset\operatorname{\mathbb{R}}^{N}$ with independent, absolutely continuous components, i.e. \xi(\omega)=\big{(}\xi_{1}(\omega),\ldots,\xi_{N}(\omega)\big{)}:\Omega\to\Xi\subset\operatorname{\mathbb{R}}^{N}. This means that for every random variable $\xi_{i}$ there exists a density function $w_{i}:\operatorname{\mathbb{R}}\to[0,\infty)$, such that $\int\limits_{\operatorname{\mathbb{R}}}w_{i}(y)~{}\mathrm{d}y=1$ and $\mathbb{P}[\xi_{i}\in A]=\int\limits_{A}w_{i}(y)~{}\mathrm{d}y$, for any $A\in\mathcal{B}(\operatorname{\mathbb{R}})$, for all $i=1,\ldots,N$. Here $\mathcal{B}(\operatorname{\mathbb{R}})$ is the Borel $\sigma$-algebra on $\operatorname{\mathbb{R}}$. Moreover, the joint density function $w$ of the random vector $\xi=(\xi_{1},\ldots,\xi_{N})$ can be written as

[TABLE]

The random vector induces a probability measure $\tilde{\mathbb{P}}(B):=\mathbb{P}(\xi^{-1}(B))$ for all $B\in\mathcal{B}(\Xi)$ on the measurable space $(\Xi,\mathcal{B}(\Xi))$. This measure is called the law of $\xi$ and in the following we work on the image probability space $(\Xi,\mathcal{B}(\Xi),\tilde{\mathbb{P}})$.

For a Banach space $E$ and its Borel $\sigma$-algebra $\mathcal{B}(E)$, we consider the weighted $L_{w}^{p}$-spaces equipped with the norms

[TABLE]

2.2 Statement of the Problem

We start with the following one-dimensional hyperbolic system of $m\in\operatorname{\mathbb{N}}$ (deterministic) conservation laws.

[TABLE]

Here, $u(t,x)\in\operatorname{\mathcal{U}}\subset\operatorname{\mathbb{R}}^{m}$ is the vector of conserved quantities, $F\in C^{2}(\operatorname{\mathcal{U}};\operatorname{\mathbb{R}}^{m})$ is the flux function, and $\mathcal{U}\subset\operatorname{\mathbb{R}}^{m}$ is the state space, which is assumed to be an open set and $T\in(0,\infty)$. We make the following assumptions on the initial condition and flux function.

1. (D1)

The initial condition satisfies $u^{0}\in{L^{1}(\operatorname{\mathbb{R}};\operatorname{\mathcal{U}})}.$ 2. (D2)

The Jacobian $\operatorname{D}F$ has $m$ distinct real eigenvalues, with each characteristic field being either genuinely nonlinear or linearly degenerate [12, Def. 3.1.1, Def. 7.5.1].

For our probabilistic equation of interest we admit in (2.2) random variations in the flux and initial datum, leading to the following one-dimensional system of $m\in\operatorname{\mathbb{N}}$ random conservation laws.

[TABLE]

Here, $u(t,x,y)\in\operatorname{\mathcal{U}}\subset\operatorname{\mathbb{R}}^{m}$ is the random vector of conserved quantities and $F(\cdot,y)\in C^{2}(\operatorname{\mathcal{U}};\operatorname{\mathbb{R}}^{m})$, $\tilde{\mathbb{P}}\text{-a.s. }y\in\Xi$ is the random flux function.

For the sake of simplicity we keep the same notation for the solution and for the flux as in (2.2) 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).

1. (R1)

The uncertain initial condition satisfies $u^{0}\in L_{w}^{1}(\Xi;L^{1}(\operatorname{\mathbb{R}};\operatorname{\mathcal{U}})).$ 2. (R2)

For almost every realization $y\in\Xi$, the Jacobian $\operatorname{D}F(\cdot,y)$ has $m$ distinct real eigenvalues, and each characteristic field is either linearly degenerate or genuinely nonlinear. Moreover, we assume that $F\in L_{w}^{2}(\Xi;C^{2}(\operatorname{\mathcal{U}};\operatorname{\mathbb{R}}^{m}))$.

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 (2.2) and establish its existence and uniqueness based on the results for the deterministic hyperbolic conservation law (2.2).

3.1 Deterministic Hyperbolic Conservation Laws

Let us consider the deterministic initial value problem (2.2). We say that a strictly convex function $\eta\in C^{2}(\operatorname{\mathcal{U}};\operatorname{\mathbb{R}})$ and a function $q\in C^{2}(\operatorname{\mathcal{U}};\operatorname{\mathbb{R}})$ form an entropy/entropy-flux pair, if they satisfy $\operatorname{D}\eta\operatorname{D}F=\operatorname{D}q$. We assume that the system (2.2) 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\in L^{1}((0,T)\times\operatorname{\mathbb{R}};\operatorname{\mathcal{U}})$ is called an entropy admissible solution of (2.2), if it satisfies the following conditions:

It is a weak solution, i.e.

[TABLE]

for all $\phi\in C_{c}^{\infty}([0,T)\times\operatorname{\mathbb{R}};\operatorname{\mathbb{R}}^{m})$. 2. 2.

It satisfies the weak entropy inequality:

[TABLE]

for all $\Phi\in C_{c}^{\infty}([0,T)\times\operatorname{\mathbb{R}};\operatorname{\mathbb{R}}_{+})$.

Theorem 3.2** ([6], Theorem 2).**

Provided (D2) holds, there exists a non-empty closed domain $\operatorname{\mathcal{D}}\subset L^{1}(\operatorname{\mathbb{R}};\operatorname{\mathcal{U}})$ of integrable functions with small total variation and a semi-group $\operatorname{\mathcal{S}}(t):[0,\infty)\times\operatorname{\mathcal{D}}\to\operatorname{\mathcal{D}}$, called Standard Riemann Semigroup (SRS), that is unique (up to its domain) and which has in particular the following properties:

1. (i)

There exists a constant $L>0$, such that

[TABLE]

for all $\overline{u},\overline{v}\in\operatorname{\mathcal{D}}$ and for all $s,t\geq 0$. 2. (ii)

For $\overline{u}\in\operatorname{\mathcal{D}}$ the function $u(t,x):=(\operatorname{\mathcal{S}}(t)\overline{u})(x)$ is an entropy admissible solution of (2.2). 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 (2.2) 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.

Theorem 3.5** ([4], Corollary 2.5).**

Let the flux function $F$ satisfy (D2) and assume

[TABLE]

for some suitable positive $M\in\operatorname{\mathbb{R}}$ and some compact set $\operatorname{\mathcal{C}}\subset\operatorname{\mathbb{R}}^{m}$. For $t>0$ we denote by $\operatorname{\mathcal{S}}(t,F)$ 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 $\tilde{F}$, satisfying (D2) and $\operatorname{\mathcal{D}}(\tilde{F})\subseteq\operatorname{\mathcal{D}}(F)$, it holds that

[TABLE]

for all $u\in\operatorname{\mathcal{D}}(\tilde{F})$.

3.2 Existence and Uniqueness of Random Entropy Solutions

We now consider the random hyperbolic conservation law (2.2) and introduce the notion of a random entropy admissible solution for (2.2). We say that $\eta\in L_{w}^{1}(\Xi;C^{2}(\operatorname{\mathcal{U}};\operatorname{\mathbb{R}}^{m}))$,  $q\in L_{w}^{1}(\Xi;C^{2}(\operatorname{\mathcal{U}};\operatorname{\mathbb{R}}^{m}))$ form a random entropy/entropy-flux pair if $\eta(\cdot,y)$ is strictly convex $\tilde{\mathbb{P}}\text{-a.s. }y\in\Xi$ and if $\eta$ and $q$ satisfy $\operatorname{D}\eta(\cdot,y)\operatorname{D}F(\cdot,y)=\operatorname{D}q(\cdot,y)$, $\tilde{\mathbb{P}}\text{-a.s. }y\in\Xi$. We assume that the random conservation law (2.2) 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 (2.2). 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\in L_{w}^{1}(\Xi;L^{1}((0,T)\times\operatorname{\mathbb{R}};\operatorname{\mathcal{U}}))$ is called a random entropy admissible solution of (2.2), if it satisfies the following conditions:

It is a weak solution:

[TABLE]

$\tilde{\mathbb{P}}\text{-a.s. }y\in\Xi$* and for all $\phi\in C_{c}^{\infty}([0,T)\times\operatorname{\mathbb{R}};\operatorname{\mathbb{R}}^{m})$.* 2. 2.

It satisfies the weak entropy inequality:

[TABLE]

$\tilde{\mathbb{P}}\text{-a.s. }y\in\Xi$* and for all $\Phi\in C_{c}^{\infty}([0,T)\times\operatorname{\mathbb{R}};\operatorname{\mathbb{R}}_{+})$.*

Remark 3.7**.**

1. (i)

Let $u$ be a random entropy admissible solution, then for almost any fixed realization $\tilde{y}\in\Xi$, the function $u(\cdot,\cdot,\tilde{y})$ is an entropy admissible solution of the deterministic version of (2.2) in the sense of Section 3.1. 2. (ii)

The function u(t,x,y):=\big{(}\operatorname{\mathcal{S}}(t,F(\cdot,y))u^{0}(\cdot,y)\big{)}(x), where $\{\operatorname{\mathcal{S}}(t,F(\cdot,y))\}_{t\geq 0}$ is the SRS from Theorem 3.2 associated with the flux-function $F(\cdot,y)$, obviously satisfies (1) and (2) $\tilde{\mathbb{P}}\text{-a.s. }y\in\Xi$. 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 (2.2).

To ensure the existence of a random entropy admissible solution of (2.2) by applying Theorem 3.2 and Theorem 3.5 path-wise in $\Xi$ we make the following assumptions:

1. (R3)

We define $\operatorname{\mathcal{D}}:=\bigcap_{y\in\Xi}\operatorname{\mathcal{D}}(F(\cdot,y))$, where $\operatorname{\mathcal{D}}(F(\cdot,y))$ is the domain of the SRS from (3.3) in Theorem 3.5. We assume that $\operatorname{\mathcal{D}}\neq\text{\O}$ and $u^{0}(\cdot,y)\in\operatorname{\mathcal{D}}$, $\tilde{\mathbb{P}}\text{-a.s. }y\in\Xi$. 2. (R4)

There exists a compact and convex set $\operatorname{\mathcal{C}}\subset\operatorname{\mathcal{U}}$, s.t. $\operatorname{\mathcal{S}}(t,F(\cdot,y))u^{0}(\cdot,y)(x)\in\operatorname{\mathcal{C}}$, a.e. $(t,x,y)\in(0,T)\times\operatorname{\mathbb{R}}\times~{}\Xi$ and $u^{0}(x,y)\in\operatorname{\mathcal{C}}$, a.e. $(x,y)\in\operatorname{\mathbb{R}}\times~{}\Xi$.

Theorem 3.8**.**

*Let the assumptions (R1)-(R4) hold. For $\tilde{\mathbb{P}}\text{-a.s. }y\in\Xi$ we define $u(t,x,y):=\operatorname{\mathcal{S}}(t,F(\cdot,y))u^{0}(\cdot,y)(x)$, where $\{\operatorname{\mathcal{S}}(t,F(\cdot,y))\}_{t\geq 0}$ is the SRS from Theorem 3.2 associated with the flux-function $F(\cdot,y)$.

Then $u$ is a random entropy admissible solution of (2.2). 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 (2.2).*

Proof.

The function $u$ is path-wise the unique SRS-induced entropy solution of (2.2) by construction. Note that we have assumed $u^{0}(\cdot,y)\in\operatorname{\mathcal{D}}\subset\operatorname{\mathcal{D}}(\cdot,F(\cdot,y))$, $\tilde{\mathbb{P}}\text{-a.s. }y\in\Xi$. It remains to show, that $u$ is a random variable, i.e.

[TABLE]

is a measurable map. To this end we define the vector space

[TABLE]

equipped with the norm

[TABLE]

Using Theorem 3.2 (i) and Theorem 3.5, which we can apply due to assumptions (R3) and (R4), we deduce

[TABLE]

for $\tilde{\mathbb{P}}$-a.s. $y,\tilde{y}\in\Xi$.

Hence, the mapping $\mathrm{S}(t):(\overline{u},\overline{}F)\ni E_{1}\to L^{1}(\operatorname{\mathbb{R}};\operatorname{\mathbb{R}}^{m})$, $\mathrm{S}(t)(\overline{u},\overline{}F):=\operatorname{\mathcal{S}}(t,\overline{F})\overline{u}(\cdot)$ is continuous for all $t>0$. Due to the finite time horizon we immediately deduce that the mapping

[TABLE]

is also continuous. Finally, it follows from assumptions (R1) and (R2) that the mapping

[TABLE]

is measurable. Thus, $u(\cdot,\cdot,y)=\operatorname{\mathcal{S}}(\cdot,F(\cdot,y))u^{0}(\cdot,y)=\mathrm{S}\circ\mathrm{S}_{0}(y)$ is a composition of measurable mappings and hence measurable itself. ∎

4 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 (2.2). 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 $\Xi$ as for example in [1, 36]. Additionally, we recall the Multi-Element method which decomposes the random space $\Xi$ 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.

4.1 The Stochastic Collocation Method

The idea of the SC method is to approximate the random entropy admissible solution of (2.2) by a polynomial interpolant in the random space, where the interpolant is supposed to satisfy (2.2) at collocation points $\{y_{i}\}_{i=0}^{M}\subset\operatorname{\mathbb{R}}$, $M\in\operatorname{\mathbb{N}}$. The exact solution $u(\cdot,\cdot,y_{i})$ at a given collocation point $y_{i}$, $i=0,\ldots,M$, is then approximated by a discrete solution $u_{h}(\cdot,\cdot,y_{i})$ using any suitable numerical method.

For a multi-dimensional random space $\Xi\subset\operatorname{\mathbb{R}}^{N}$, we define $\Xi_{i}:=\xi_{i}(\Omega)$ and follow [2] to first define the space $\mathcal{P}_{q}(\Xi)$ of tensor product polynomials of maximal degree $q\in\operatorname{\mathbb{N}}_{0}$ by

[TABLE]

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\in\operatorname{\mathbb{N}}_{0}$ for every random dimension we could also consider variable polynomial degrees $q_{i}\in\operatorname{\mathbb{N}}$, $i=1,\ldots,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 $\mathcal{K}:=\{\mathrm{k}=(k_{1},\ldots,k_{N})^{\top}\in\operatorname{\mathbb{N}}_{0}^{N}:k_{i}\leq q,~{}i=1\ldots,N\}$ be the corresponding multi-index set and define the collocation points $y_{\mathrm{k}}=(y_{1,k_{1}},\ldots,y_{N,k_{N}})\in\Xi$, $\mathrm{k}\in\operatorname{\mathcal{K}}$. As a basis of $\mathcal{P}_{q}(\Xi_{i})$ we choose the Lagrange basis $\{l_{i,k}\}_{k=0}^{q}$ associated with the collocation points $\{y_{i,k}\}_{k=0}^{q}$, such that

[TABLE]

for all $i=1,\ldots,N$. We then define the multivariate Lagrange polynomials via

[TABLE]

Using the collocation points $\{y_{\mathrm{k}}\}_{\mathrm{k}\in\operatorname{\mathcal{K}}}$ as input parameters in (2.2) yields $\text{card}(\operatorname{\mathcal{K}})=(q+1)^{N}$ (uncoupled) collocated initial value problems:

[TABLE]

Here and in the following we consider $\Lambda\in\{[0,1]_{per},\operatorname{\mathbb{R}}\}$.

Remark 4.2**.**

The well-posedness result in Theorem 3.8 only covers $\Lambda=\operatorname{\mathbb{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 $\Lambda=[0,1]_{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_{\mathrm{k}}$ we denote the corresponding numerical approximation at time $t_{n}=t_{n}(y_{\mathrm{k}})$ by $u_{h}^{n}(\cdot,y_{\mathrm{k}}):=u_{h}(t_{n},\cdot,y_{\mathrm{k}})$. Let us assume that the time-partition $\{t_{n}\}_{n=0}^{N_{t}}$ is the same for every collocation point $\{y_{\mathrm{k}}\}_{\mathrm{k}\in\operatorname{\mathcal{K}}}$. The numerical approximation of (2.2) at time $t=t_{n}$ can then be written as

[TABLE]

An important aspect of the SC method is the choice the collocation points $\{y_{i,k}\}_{k=0}^{q}\subset\Xi_{i}$. Depending on the distribution of the random variable $\xi_{i}$ we choose the collocation points as zeros of the corresponding (orthogonal) chaos polynomials [37]. For example, if $\xi_{i}\sim\mathcal{U}(a,b)$ is uniformly distributed, we choose $\{y_{i,k}\}_{k=0}^{q}$ to be the roots of the $(q+1)$-th Legendre polynomial. For a Gaussian distribution we use the roots of the Hermite polynomials accordingly.

One can then approximate the mean of $u_{h}^{n}(x,\cdot)$ via numerical quadrature, i.e.

[TABLE]

Here $w_{\mathrm{k}}$ are products of the corresponding one-dimensional weights.

4.2 The Multi-Element Stochastic Collocation Method

A major drawback of any global approximation approach in $\Xi$ for hyperbolic conservation laws is that, due to the Gibbs phenomenon, the interpolant may oscillate for discontinuous solutions, cf. [33, 28]. To overcome this issue, we employ the Multi-Element (ME) approach as presented in [34], i.e. we subdivide $\Xi$ into disjoint elements and consider a local approximation of (2.2) on every random element.

For the ease of presentation we assume that $\Xi=[0,1]^{N},$ and let $0=d_{1}<d_{2}<\ldots<d_{N_{\Xi}+1}=1$ be a decomposition of $[0,1]$. We define $D_{n}:=[d_{n},d_{n+1})$, for $n=1,\ldots,N_{\Xi}-1,$ and $D_{N_{\Xi}}:=[d_{N_{\Xi}},d_{N_{\Xi}+1}]$. Introducing the set $\mathcal{M}:=\{\mathrm{m}=(m_{1},\ldots,m_{N})^{\top}\in\operatorname{\mathbb{N}}_{0}^{N}:m_{i}\leq N_{\Xi},~{}i=1\ldots,N\}$ allows us to define for $\mathrm{m}\in\mathcal{M}$, the Multi-Element D_{\mathrm{m}}:=\bigtimes\limits_{i=1}^{N}D_{m_{i}}. Hence, we consider a new local random variable $\xi^{\mathrm{m}}:\xi^{-1}(D_{\mathrm{m}})\to D_{\mathrm{m}}$ on the local probability space $(\xi^{-1}(D_{\mathrm{m}}),\operatorname{\mathcal{F}}\cap\xi^{-1}(D_{\mathrm{m}}),\mathbb{P}(\cdot|\xi^{-1}(D_{\mathrm{m}})))$. Using Bayes’ rule we can compute the local probability density function of $\xi^{\mathrm{m}}$ via

[TABLE]

where $\mathbb{P}(\xi^{-1}(D_{\mathrm{m}}))>0$ for $\mathrm{m}\in\mathcal{M}$ can be assumed w.l.o.g., due to the independence of the corresponding random variables. We parametrize the uncertain input in (2.2) using the local random variable $\xi^{\mathrm{m}}$ and consider a local approximation on every $D_{\mathrm{m}}$ at time $t=t_{n}$, $n=0,\ldots,N_{t}$,

[TABLE]

for all $\mathrm{m}\in\mathcal{M}$. Here, $\{y_{\mathrm{k}}^{\mathrm{m}}\}_{\mathrm{k}\in\operatorname{\mathcal{K}}}\subset D_{\mathrm{m}}$ are the local collocation points in $D_{\mathrm{m}}$, $\mathrm{m}\in\mathcal{M}$. The global Multi-Element Stochastic Collocation (ME-SC) approximation at time $t=t_{n}$, can then be written as

[TABLE]

where $\chi_{D_{\mathrm{m}}}$ is the indicator function of $D_{\mathrm{m}}$.

4.3 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_{\mathrm{k}}$, $\mathrm{k}\in\operatorname{\mathcal{K}}$, a computable space-time reconstruction $\hat{u}^{st}(y_{\mathrm{k}})=\hat{u}^{st}(\cdot,\cdot,y_{\mathrm{k}})\in W_{\infty}^{1}((0,T);V_{p+1}^{s}\cap C^{0}(\Lambda))$ of the numerical approximation $u_{h}(\cdot,\cdot,y_{\mathrm{k}})$, where $V_{p+1}^{s}$ denotes the space of piece-wise polynomials of degree $p+1$ on a triangulation of $\Lambda$. This allows us to define the space-time residual as follows.

Definition 4.3** (Space-time residual).**

We call the function $R^{st}(y_{\mathrm{k}}):=R^{st}(\cdot,\cdot,y_{\mathrm{k}})\in L_{2}((0,T)\times\Lambda;\operatorname{\mathbb{R}}^{m})$, defined by

[TABLE]

the space-time residual associated with the collocation point $y_{\mathrm{k}}$, for all $\mathrm{k}\in\operatorname{\mathcal{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 \hat{u}^{sts}\in\mathcal{P}_{q}(\Xi)\otimes\big{(}W_{\infty}^{1}(0,T);V_{p+1}^{s}\cap C^{0}(\Lambda)\big{)} defined by

[TABLE]

the space-time-stochastic reconstruction of the numerical approximation $u_{h}$ of (2.2) (see (4.1)).

Similar to the space-time reconstruction, we may plug $\hat{u}^{sts}$ into the random conservation law (2.2) 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}\in L_{w}^{2}(\Xi;L^{2}((0,T)\times\Lambda;\operatorname{\mathbb{R}}^{m}))$ by

[TABLE]

We also need this residual for the upcoming error analysis.

5 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 $\operatorname{\mathcal{C}}$, there exist $\tilde{\mathbb{P}}\text{-a.s. }y\in\Xi$ constants $0<C_{\overline{F}}(y)<\infty$ and $0<C_{\underline{\eta}}(y)<C_{\overline{\eta}}(y)<\infty$, such that,

[TABLE]

Here, for a generic function $f$, $H_{u}f$ denotes its Hessian matrix which contains all second order derivatives with respect to $u$.

5.1 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 (2.2). Let the reconstruction $\hat{u}^{sts}$ only take values in $\mathcal{C}$. Then, the difference between $u$ and the numerical solution $u_{h}^{n}$ from (4.1) satisfies

[TABLE]

for all $n=0,\ldots,N_{t}$ and for any $\tilde{\mathbb{P}}$-measurable set $\tilde{\Xi}\subseteq\Xi$. Here

[TABLE]

Proof.

Thanks to the path-wise structure we can apply the deterministic setting from [15, Theorem 5.5 ] for almost any $y\in\tilde{\Xi}\subseteq\Xi$. Integrating over $\tilde{\Xi}$ gives the desired estimate directly. ∎

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 $\mathcal{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

[TABLE]

with

[TABLE]

$R^{det}$* and $R^{stoch}$ are called the deterministic and stochastic residual.*

Proof.

For every collocation point $y_{\mathrm{k}},$ $\mathrm{k}\in\operatorname{\mathcal{K}}$, we compute in Appendix A the space-time reconstruction $\hat{u}^{st}(\cdot,\cdot,y_{\mathrm{k}})$ which fulfills

[TABLE]

Moreover, we know from (4.6) that the space-time-stochastic reconstruction $\hat{u}^{sts}=\sum\limits_{\mathrm{k}\in\operatorname{\mathcal{K}}}\hat{u}^{st}(t,x,y_{\mathrm{k}})l_{\mathrm{k}}(y)$ satisfies the relation

[TABLE]

Multiplying (5.6) by $l_{\mathrm{k}}$ and summing over $\mathrm{k}\in\operatorname{\mathcal{K}}$ yields

[TABLE]

Inserting (5.8) into (5.7) yields

[TABLE]

∎

To simplify Theorem 5.1 let us assume that the eigenvalues of the Hessian $H_{u}\eta(u,y)$ are bounded from above and below by positive numbers, for any $u\in\operatorname{\mathcal{C}}$ uniformly in $\Xi$. We let $C_{\overline{\eta}}:=\operatorname{ess}\sup\limits_{y\in\Xi}C_{\overline{\eta}}(y)<\infty$, $C_{\underline{\eta}}:=\operatorname{ess}\inf\limits_{y\in\Xi}C_{\overline{\eta}}(y)>0$ and $C_{\overline{F}}:=\operatorname{ess}\sup\limits_{y\in\Xi}C_{\overline{F}}(y)<\infty$. In our numerical examples, where we consider the compressible Euler equations, the dependence of the flux function $F$ and $\eta$ on $y$ is explicitly known and we can compute the constants $C_{\underline{\eta}},C_{\overline{\eta}},C_{\overline{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 (2.2). Then, the difference between $u$ and the numerical solution $u_{h}^{n}$ from (4.1) satisfies

[TABLE]

for $n=0,\ldots,N_{t}$ and for all $\tilde{\mathbb{P}}$-measurable sets $\tilde{\Xi}\subseteq\Xi$. Here,

[TABLE]

Remark 5.4**.**

1. (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 space-time 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. 2. (ii)

*In order for the upcoming space-stochastic adaptive algorithm based on $\mathcal{E}^{det}$, $\mathcal{E}^{stoch}$ to be efficient, we need $\mathcal{E}^{det}$ to depend solely on the space-time discretization and to be independent of the stochastic discretization. Similarly, we need $\mathcal{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 $\mathcal{E}^{det}$ is indeed unaffected by stochastic refinement. Our numerical experiment in Section 6.3.1 also shows that $\mathcal{E}^{stoch}$ is unaffected by spatial refinement. 3. (iii)

The scaling properties of $\mathcal{E}^{det}$ , resp. $R^{st}(y_{\mathrm{k}})$, were studied in **[13]**. Currently we are not able to prove any of the scaling properties of $\mathcal{E}^{stoch}$ w.r.t. to $q$ and the number of Multi-Elements. However, our numerical experiments show that $\mathcal{E}^{stoch}$ scales as desired, i.e. $\mathcal{E}^{stoch}$ shows the same qualitative behavior as the stochastic interpolation error of the exact solution. 4. (iv)

As described in **[13]** and **[16]**, $R^{det}$ scales with $\frac{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 (2.2). However, as $R^{det}$ precisely captures the positions of rarefaction waves, contact discontinuities and shocks we use $R^{det}$ and $R^{stoch}$, resp. $\mathcal{E}^{det}$ and $\mathcal{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 $\Xi$).**

As noted above, the collocation points $y_{\mathrm{k}}$ are chosen to be the zeros of the corresponding orthogonal polynomial depending on the distribution of $\xi$. The deterministic residual $R^{det}$ from (5.4) consists of Lagrange polynomials associated with the corresponding collocation points, thus Gaussian quadrature in $\Xi$ yields

[TABLE]

Hence, $\mathcal{E}^{det}$ inherits the convergence order of $R^{st}(y_{\mathrm{k}})$.

5.2 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_{\mathrm{k}}$, $\mathrm{k}\in\operatorname{\mathcal{K}}$. We slightly abuse the notation from (5.9) and write for every physical cell $I\in\mathcal{T}_{n}$, $\mathcal{E}^{det}_{\mathrm{k}}(t_{n},t_{n+1},I):=\|R^{st}(y_{\mathrm{k}})\|_{L^{2}((t_{n},t_{n+1})\times I)}$, which is the cell-wise indicator for the spatial refinement in $\Lambda$.

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 $\mathcal{E}^{det}_{\mathrm{k}}(t_{n},t_{n+1},I)$ on every cell $I\in\mathcal{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 $\mathcal{E}^{det}_{\mathrm{k}}(t_{n},t_{n+1},I)$, 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\in\operatorname{\mathbb{N}}$, to restrict the fineness of the adaptive mesh. The algorithm reads as follows:

Remark 5.6**.**

After every projection step in line three of Algorithm 1 we apply the TVBM slope limiter $\Lambda\Pi_{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 $\mathcal{E}^{stoch}$ as local indicator for stochastic refinement.

6 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 $\mathcal{E}^{det}$ and $\mathcal{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 RK-method 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

[TABLE]

or the Lax-Friedrichs numerical flux

[TABLE]

In our example, the uncertainty is uniformly distributed. Therefore, we use the zeros of the Gauß-Legendre polynomials as collocation points. Computing $\mathcal{E}^{det},\mathcal{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 (2.2) the one-dimensional compressible Euler equations for the flow of an ideal gas, which are given by

[TABLE]

where $\rho$ describes the mass density, $m$ the momentum and $E$ the energy of the gas. The constitutive law for pressure $p$ reads

[TABLE]

with the adiabatic constant $\gamma=1.4$ if not specified otherwise. In the following figures we refer to the quantity $\|m(T,\cdot,\cdot)-m_{h}^{N_{t}}(\cdot,\cdot)\|_{L_{w}^{2}(\Xi;L^{2}(\Lambda))}$ at final computational time $T$ as numerical error, unless otherwise stated. We also plot the residuals $\mathcal{E}^{det}(T)$ and $\mathcal{E}^{stoch}(T)$ 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),

[TABLE]

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.

6.1 A Deterministic Problem with Smooth Solution

In this numerical example, we study the scaling properties of the deterministic residual $\mathcal{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

[TABLE]

The numerical solution is computed on the domain $\Lambda=[0,1]_{\mathrm{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 $\mathcal{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.

6.2 Deterministic Adaptivity: Sod Shock Tube Problem

In this numerical experiment we apply the adaptive spatial mesh refinement from Algorithm 1 to the Sod shock tube problem. The Riemann data for this problem is given by

[TABLE]

The numerical solution is computed on the domain $\Lambda=[0,1]$ 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 (1) we compare the $L^{1}(\Lambda)$- and $L^{2}(\Lambda)$-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 (1).

Remark 6.2**.**

As discussed in Remark 5.4 (iv), $R^{det}$ scales with $\frac{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].

6.3 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 $\Xi$ and a random flux function.

6.3.1 A One-Dimensional Random Space, $q$-Refinement

We modify the exact solution from Section 6.1 in the following way,

[TABLE]

The numerical solution is computed on $\Lambda=[0,1]_{\mathrm{per}}$ up to $T=0.2$, the uncertainty $y$ stems from an uniform distribution, i.e. $\xi\sim\mathcal{U}(0,8)$. 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 $\Xi$, 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 $\mathcal{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 $\mathcal{E}^{stoch}$ is not altered by a finer or coarser space discretization. Moreover, the deterministic residual $\mathcal{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.

6.3.2 Mesh Refinement in $\Xi$ and Random Flux Function

In this example we examine the scaling properties of $\mathcal{E}^{stoch}$ under mesh refinements for a two-dimensional random space $\Xi\subset\operatorname{\mathbb{R}}^{2}$. We consider the same smooth function as in Section 6.3.1,

[TABLE]

with $\xi_{1}\sim\mathcal{U}(0,8)$. Moreover, we consider a random adiabatic constant. We assume that $\gamma=\xi_{2}\sim\mathcal{U}(1.4,1.6)$ 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 $\Lambda=[0,1]_{\mathrm{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\in\{1,2\}$.

Figure (3) illustrates the behavior of the stochastic residual $\mathcal{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\in\{1,2\}$, the stochastic residual converges approximately with the expected rate of convergence, which is $(q+1)/2$, cf. [34]. Like for the $q$-refinement in Section 6.3.1, the deterministic residual $\mathcal{E}^{det}$ stays constant, when we increase the number of MEs.

6.4 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,

[TABLE]

where

[TABLE]

is a discontinuous amplitude. For the spatial domain $\Lambda=[0,1]_{\mathrm{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 $\xi_{1},\xi_{2}\sim\operatorname{\mathcal{U}}(0,1)$. For the ME-SC method we consider a linear interpolant, i.e. $q=1$.

In Figure (4) we plot the error and the spatial resp. stochastic residual versus the number of MEs and in Figure (4) we show the error of the uniform and adaptive method versus cpu time. In Figure (4) we can observe that for the uniform stochastic refinement, both the error and the stochastic residual $\mathcal{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 (4), where we reach an error reduction in significantly less time compared to uniform refinement.

6.5 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

[TABLE]

where $\xi_{1},\xi_{2}\sim\mathcal{U}(-1,1)$ and $\Lambda=[0,1]$. 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_{w}^{2}(\Xi;L^{2}(\Lambda))}$-norm. Note that we do not have an exact solution at hand for this problem, but due to Jensen’s inequality,

[TABLE]

The reference expectation $\operatorname{\mathbb{E}}(u(T,\cdot,\cdot))$ is computed using a Monte-Carlo method with an exact Riemann solver with 500000 samples.

In Figure (5) we show the numerical error as in (6.7) and we also consider the error $\|\operatorname{\mathbb{E}}(u(T,\cdot,\cdot))-\operatorname{\mathbb{E}}(u_{h}^{N_{t}}(\cdot,\cdot))\|_{L^{1}(\Lambda)}$ for an increasing number of MEs, i.e. for increasing $N_{\Xi}$. 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 (5)), 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.

Appendix A RKDG Method and Space-Time Reconstructions

In the following we describe the space-time discretization of (CIVP) that we use and the space-time reconstruction of the numerical solution.

A.1 The RKDG Method

For the space-time discretization of (CIVP) we use a Runge–Kutta DG method (RKDG). We recall the DG spatial discretization as for example in [10]. For the ease of presentation we neglect the dependence of the flux $F$, the spatial mesh and the DG spaces on the collocation points $\{y_{\mathrm{k}}\}_{\mathrm{k}\in\operatorname{\mathcal{K}}}$.

Let $\mathcal{T}:=\{I_{k}\}_{k=0}^{{N_{s}-1}}$, $I_{k}:=(x_{k},x_{k+1})$ be a quasi-uniform triangulation of $\Lambda=[0,1]_{\mathrm{per}}$. We set $h_{k}=(x_{k+1}-x_{k})$, $h_{\max}=\max\limits_{k}h_{k}$, $h_{\min}=\min\limits_{k}h_{k}$ for the spatial mesh and identify $x_{0}=x_{N_{s}}$ to account for periodic boundary conditions. Further let $0=t_{0}<t_{1}<\ldots<t_{N_{t}}=T$ be a temporal decomposition of $[0,T]$ and define $\Delta t_{n}:=(t_{n+1}-t_{n})$, $\Delta t=\max\limits_{n}~{}\Delta t_{n}$. With each time-interval $(t_{n},t_{n+1}]$ we associate a (possibly different) partition $\mathcal{T}_{n}$ and associated DG space

[TABLE]

With $\mathcal{L}_{V_{p,n}^{s}}$ we denote the $L^{2}$-projection mapping into the DG space $V_{p,n}^{s}$.

Following [14] we call the function $u_{h}$ a generalized semi-discrete DG approximation of (CIVP) if it satisfies for $u_{h}^{-1}:=\mathcal{L}_{V_{p,0}^{s}}u^{0}$ the following equations. For every $n=0,\ldots,N_{t}$, $u_{h}^{n}\mid_{[t_{n},t_{n+1}]}\in C^{1}((t_{n},t_{n+1});V_{p,n}^{s})\cap C^{0}([t_{n},t_{n+1}];V_{p,n}^{s})$,

[TABLE]

where $L_{h}^{n}:V_{p,n}^{s}\to V_{p,n}^{s}$ is defined by

[TABLE]

The numerical solution $u_{h}$ is defined through $u_{h}(0):=u_{h}^{-1}$ and $u_{h}\mid_{(t_{n},t_{n+1}]}:=u_{h}^{n}\mid_{(t_{n},t_{n+1}]}$.

Here, $\operatorname{G}:\operatorname{\mathbb{R}}^{m}\times\operatorname{\mathbb{R}}^{m}\to\operatorname{\mathbb{R}}^{m}$ denotes a numerical flux, the spatial traces are defined as $\psi(x^{\pm}):=\lim\limits_{h\searrow 0}\psi(x\pm h)$ and $\operatorname{[\![}\psi_{h}\operatorname{]\!]}_{i}:=(\psi_{h}(x_{i}^{-})-\psi_{h}(x_{i}^{+}))$ are jumps.

The initial-value problem (DG) can now be solved numerically by any single- or multi-step method. We focus on $K$-th order Runge-Kutta time-step methods as in [11, 22]. Furthermore, $\Lambda\Pi_{h}:\operatorname{\mathbb{R}}^{m}\to\operatorname{\mathbb{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_{h}^{n}(t_{n})\in V_{p,n}^{s}$ reads as follows:

The parameters $\alpha_{jl}$ satisfy the conditions $\alpha_{jl}\geq 0$, $\sum\limits_{l=0}^{j-1}\alpha_{jl}=1$ , and if $\beta_{jl}\neq 0$, then $\alpha_{jl}\neq 0$ for all $j=1,\ldots,S$, $l=0,\ldots,j$.

A.2 Space-Time Reconstruction

Our analysis relies on reconstructing the numerical solution $\{u_{h}^{n}\}_{n=0}^{N_{t}}$ to a Lipschitz continuous function in space and time. We structure the reconstruction process as follows:

Computation of a temporal reconstruction $\hat{u}^{t}$:

We first compute the temporal reconstruction as proposed in [13].

Let $\{u_{h}^{0},\ldots,u_{h}^{N_{t}}\}$ be a sequence of approximate solutions of (DG) at points $\{t_{n}\}_{n=0}^{N_{t}}$ in time, where we assume that all approximate solutions are interpolated onto a reference mesh $\mathcal{T}$, which is a common refinement of all meshes. With $V_{p}^{s}$ we denote the DG space associated with $\mathcal{T}$, hence $u_{h}(t_{n})\in V_{p}^{s}$ for all $n=0,\ldots,N_{t}$.

For the reconstruction in time we define the spaces of piecewise polynomials in time of degree $r$ by

[TABLE]

Using Hermite interpolation on each time interval $[t_{n},t_{n+1}]$, we construct the temporal reconstruction $\hat{u}^{t}\in V_{r}^{t}((0,T);V_{p}^{s})$. 2. 2.

Computation of a space-time reconstruction $\hat{u}^{st}$ using the time reconstruction $\hat{u}^{t}$:

With the temporal reconstruction $\hat{u}^{t}$ at hand, we define the space-time reconstruction $\hat{u}^{st}$ of the DG-solutions of (DG). The analysis in [13] requires numerical fluxes $\operatorname{G}$ which admit a special representation. In particular, there needs to exist a locally Lipschitz function $w:\operatorname{\mathcal{U}}\times\operatorname{\mathcal{U}}\to\operatorname{\mathcal{U}}$, with the additional property $w(u,u)=u$, such that $\operatorname{G}$ can either be expressed as

[TABLE]

or as

[TABLE]

where $\nu\in\operatorname{\mathbb{N}}$ and for some matrix-valued function $\mu$, which has the property that for any compact $K\subset\operatorname{\mathcal{U}}$ there exists a $\mu_{K}>0$, such that $|\mu(u,v;h)|\leq\mu_{K}(1+\frac{|v-u|}{h})$, for $h$ small enough.

Remark A.1**.**

For our numerical computations we consider the following numerical fluxes.

• •

The Lax-Wendroff flux: $\operatorname{G}(u,v)=F(w(u,v))$ with $w(u,v)=\frac{u+v}{2}-\frac{\Delta t}{2h}(F(u)-F(v))$, satisfies (A.2).

• •

The Lax-Friedrichs flux : G(u,v)=\frac{1}{2}\Big{(}F(u)+F(w)\Big{)}+\lambda(w-u) satisfies (A.3), with $\nu=0$, $w(u,v):=\frac{1}{2}(u+v)$ and $\mu(u,v;h):=\lambda I-\frac{F(u)-2F(w(u,v))+F(v)}{2|v-u|^{2}}\otimes(u-v)$.

We define the spatial reconstruction which is applied to the temporal reconstruction $\hat{u}^{t}(t,\cdot)$ for each $t\in(0,T)$ using the function $w$ (cf. [13, 15]).

Definition A.2** (Space-time reconstruction).**

Let $\hat{u}^{t}$ be the temporal reconstruction of a sequence $\{u_{h}^{n}\}_{n=0}^{N_{t}}$ of solutions of the fully discrete scheme of (DG) using a numerical flux satisfying (A.2) or (A.3). The space-time reconstruction $\hat{u}^{st}(t,\cdot)\in V_{p+1}^{s}$ is defined as the solution of

[TABLE]

We have the following property of the space-time reconstruction.

Lemma A.3** ([13], Lemma 24).**

Let $\hat{u}^{st}$ be the space-time reconstruction from Definition A.2. For each $t\in(0,T)$, the function $\hat{u}^{st}(t,\cdot)$ is well defined. Moreover,

[TABLE]

.