A Novel Weakly-Intrusive Non-linear Multiresolution Framework for Uncertainty Quantification in Hyperbolic Partial Differential Equations

Abstract

In this paper, a novel multiresolution framework, namely the Truncate and Encode (TE) approach, previously proposed by some of the authors (Abgrall et al. in J Comput Phys 257:19–56, 2014. doi:10.1016/j.jcp.2013.08.006), is generalized and extended for taking into account uncertainty in partial differential equations (PDEs). Innovative ingredients are given by an algorithm permitting to recover the multiresolution representation without requiring the fully resolved solution, the possibility to treat a whatever form of pdf and the use of high-order (even non-linear, i.e. data-dependent) reconstruction in the stochastic space. Moreover, the spatial-TE method is introduced, which is a weakly intrusive scheme for uncertainty quantification (UQ), that couples the physical and stochastic spaces by minimizing the computational cost for PDEs. The proposed scheme is particularly attractive when treating moving discontinuities (such as shock waves in compressible flows), even if they appear during the simulations as it is common in unsteady aerodynamics applications. The proposed method is very flexible since it can easily coupled with different deterministic schemes, even with high-resolution features. Flexibility and performances of the present method are demonstrated on various numerical test cases (algebraic functions and ordinary differential equations), including partial differential equations, both linear and non-linear, in presence of randomness. The efficiency of the proposed strategy for solving stochastic linear advection and Burgers equation is shown by comparison with some classical techniques for UQ, namely Monte Carlo or the non-intrusive polynomial chaos methods.

Appendix: The High-Order MUSCL Hancock Method (MHM)

In this section the MUSCL Hancock Method (MHM) is briefly recalled. A very interesting and exhaustive presentation of this method could be found in [38].

Let us consider a 1D scalar conservation law

$$\begin{aligned} \dfrac{\partial u(x,t)}{\partial t} + \dfrac{\partial f(u(x,t))}{\partial x} = 0, \end{aligned}$$

(46)

where $x \in \varOmega \subset {\fancyscript{R}}$ is the physical space and $t \in T \subset {\fancyscript{R}}^+$ is the time space. In the context of finite volume scheme the physical space is divided in a set of non-overlapping cells ${{\fancyscript{C}}}_i$ with $\varOmega = \bigcup _i {{\fancyscript{C}}}_i$. The classical first order Godunov scheme, applied to (46), is obtained introducing the so-called cell-average $\bar{u}_i$ on each cell ${{\fancyscript{C}}}_i$:

$$\begin{aligned} \bar{u}_i(t) = \frac{1}{ |{{\fancyscript{C}}}_i| } \int _{{{\fancyscript{C}}}_i}{ u(x,t) \mathrm {d}x }, \end{aligned}$$

(47)

where $|{{\fancyscript{C}}}_i|$ indicated the volume of the cell.

After the integration on each cell ${{\fancyscript{C}}}_i$, it can be written

$$\begin{aligned} |{\fancyscript{C}}_i| \dfrac{{\mathrm {d}} \bar{u}_i}{{\mathrm {d}}t} + f(u(x_L,t)) - f(u(x_R,t)) = 0, \end{aligned}$$

(48)

with ${\fancyscript{C}}_i = [ x_{i_L}, x_{i_R} ]$ and where $x_{i_L}$ and $x_{i_R}$ indicate the left and right interfaces.

Integrating in time the Eq. (48) between $t_n$ and $t_{n+1} = t_n +\varDelta t$, it follows that

$$\begin{aligned}&|{\fancyscript{C}}_i| \left( \bar{u}_i^{n+1} - \bar{u}_i^{n} \right) + \int _{t_n}^{t_{n+1}}{f(u(x_L,t)) \mathrm {d}t} - \int _{t_n}^{t_{n+1}}{f(u(x_R,t)) \mathrm {d}t} \nonumber \\&\quad =|{\fancyscript{C}}_i| \left( \bar{u}_i^{n+1} - \bar{u}_i^{n} \right) + \varDelta t \left( F_{i_L}^n - F_{i_R}^n \right) = 0, \end{aligned}$$

(49)

where a numerical approximation for the flux along the interface $x_{i_L}$ (and $x_{i_R}$) holds

$$\begin{aligned} F_{i_L}^n \approx \frac{1}{\varDelta t} \int _{t_n}^{t_{n+1}}{f(u(x_{i_L},t)) \mathrm {d}t}. \end{aligned}$$

(50)

The final form for the first order Godunov scheme is

$$\begin{aligned} \bar{u}_i^{n+1} = \bar{u}_i^n - \frac{\varDelta t}{|{\fancyscript{C}}_i|} \left( F_{i_L}^n - F_{i_R}^n \right) . \end{aligned}$$

(51)

As pointed out by LeVeque in [34] for hyperbolic problem, the information propagates with finite speed and it is reasonable to suppose that each numerical flux, at the interface, is function of the solution on the two cells to which it belongs: $F_{i_L}^n = F_{i_L}^n(\bar{u}_{i-1}^n,\bar{u}_{i}^n)$ and $F_{i_R}^n = F_{i_R}^n(\bar{u}_{i}^n,\bar{u}_{i+1}^n)$.

In this work, an exact Riemann solver is used to compute the numerical flux. In particular, given two constant states $\bar{u}_{i-1}^n$ and $\bar{u}_{i}^n$ separated by the interface, the exact solution of the problem, the so-called Riemann problem, can be found and the solution at the interface computed. The numerical flux is then equal to the flux function evaluated by knowing the exact solution at the interface. In the following, the numerical flux function obtained via the solution of the exact Riemann problem is indicated as $F_{i_L}^n = {\fancyscript{F}}^{\mathrm {RM}}(u_{i-1}^n, u_{i}^n)$ for the left interface or $F_{i_R}^n = {\fancyscript{F}}^{\mathrm {RM}}(u_{i}^n, u_{i+1}^n)$ for the right interface.

The first order Godunov scheme introduces a great amount of numerical diffusion and then van Leer [34, 38] proposed to consider non constant data on each cell to achieve a higher accuracy. From this idea, the so-called Monotone Upstream-centred Scheme for Conservation Laws (MUSCL) has been proposed. In this work, a piecewise linear approximation is used for the solution $u(x,t)$ on the cell $|{\fancyscript{C}}_i|$:

$$\begin{aligned} u(x,t_n) = \bar{u}_i^n + \sigma _i^n ( x-x_i ) \quad \mathrm {with} \quad x_{i_L} \le x \le x_{i_R}, \end{aligned}$$

(52)

in which $\sigma _i^n$ is the slope. The choice of $\sigma _i^n=0$ lead to the Godunov scheme. Computing the slope $\sigma _i^n$ as a function of only the cell averaged solution in the neighboring cells, i.e. $\sigma _i^n = \sigma _i^n(\bar{u}_{i-1}^n, \bar{u}_{i+1}^n)$, is not the best choice. If centered slope are used, spurious oscillations can be introduced in discontinuous solution. In practice, a slope limiter should be introduced near the discontinuity to avoid oscillations. In this work the Roe’s superbee limiter is employed in which

$$\begin{aligned} \left\{ \begin{array}{ll} \sigma _i^n &{}= \mathrm {maxmod}\left( \sigma _{(1)}^n, \sigma _{(2)}^n \right) \\ \sigma _{(1)}^n &{}= \mathrm {minmod}\left( \left( \frac{\bar{u}_{i+1}^n - \bar{u}_{i}^n }{|{\fancyscript{C}}_i|}\right) , 2\left( \frac{\bar{u}_{i}^n - \bar{u}_{i-1}^n }{|{\fancyscript{C}}_i|}\right) \right) \\ \sigma _{(2)}^n &{}= \mathrm {minmod}\left( 2\left( \frac{\bar{u}_{i+1}^n - \bar{u}_{i}^n }{|{\fancyscript{C}}_i|}\right) , \left( \frac{\bar{u}_{i}^n - \bar{u}_{i-1}^n }{|{\fancyscript{C}}_i|}\right) \right) , \\ \end{array} \right. \end{aligned}$$

(53)

where the $\mathrm {minmod}$ and $\mathrm {maxmod}$ are defined as follows

$$\begin{aligned} \mathrm {minmod}(a,b)= & {} \left\{ \begin{array}{ll} a \quad \mathrm {if} \quad &{} |a|<|b| \quad \mathrm {and} \quad ab > 0 \\ b \quad \mathrm {if} \quad &{} |a|>|b| \quad \mathrm {and} \quad ab > 0 \\ 0 \quad \mathrm {if} \quad &{} ab <= 0 \end{array} \right. \nonumber \\ \mathrm {maxmod}(a,b)= & {} \left\{ \begin{array}{ll} a \quad \mathrm {if} \quad &{} |a|>|b| \quad \mathrm {and} \quad ab > 0 \\ b \quad \mathrm {if} \quad &{} |a|<|b| \quad \mathrm {and} \quad ab > 0 \\ 0 \quad \mathrm {if} \quad &{} ab <= 0. \end{array} \right. \end{aligned}$$

(54)

The fully discrete second order MHM, to compute the cell averaged solution $\bar{u}_i^{n+1}$, consists of the following three steps:

• For each cell ${\fancyscript{C}}_{\ell } \in \left\{ {\fancyscript{C}}_{i-1}, {\fancyscript{C}}_{i}, {\fancyscript{C}}_{i+1}\right\}$ the solution at the interface is computed according to

  $$\begin{aligned} \left\{ \begin{array}{ll} u_{\ell _L}^n &{}= \bar{u}_{\ell }^n - \sigma _{\ell }^n \frac{|{\fancyscript{C}}_{\ell }|}{2} \\ u_{\ell _R}^n &{}= \bar{u}_{\ell }^n + \sigma _{\ell }^n \frac{|{\fancyscript{C}}_{\ell }|}{2} \\ \end{array} \right. \end{aligned}$$

  (55)

• On each cell ${\fancyscript{C}}_\ell \in \left\{ {\fancyscript{C}}_{i-1}, {\fancyscript{C}}_{i}, {\fancyscript{C}}_{i+1} \right\}$ the solution is evolved using an half time step:

  $$\begin{aligned} \left\{ \begin{array}{ll} u_{\ell _R}^{\Uparrow } &{}= u_{\ell _R}^n + \frac{1}{2} \frac{\varDelta t}{|{\fancyscript{C}}_{\ell }|} \left( f(u_{{\ell }_L}^n) - f(u_{{\ell }_R}^n) \right) \\ u_{{\ell }_L}^{\Uparrow } &{}= u_{\ell _L}^n + \frac{1}{2} \frac{\varDelta t}{|{\fancyscript{C}}_{\ell }|} \left( f(u_{{\ell }_L}^n) - f(u_{{\ell }_R}^n) \right) \\ \end{array} \right. \end{aligned}$$

  (56)

• The cell averaged value on the cell ${\fancyscript{C}}_i$ is evolved following

  $$\begin{aligned} \bar{u}_{i}^{n+1} = \bar{u}_{i}^{n} - \frac{\varDelta t}{|{\fancyscript{C}}_i|} \left( {\fancyscript{F}}^{\mathrm {RM}}\left( u_{{i-1}_R}^{\Uparrow }, u_{i_L}^{\Uparrow } \right) - {\fancyscript{F}}^{\mathrm {RM}}\left( u_{i_R}^{\Uparrow }, u_{{i+1}_L}^{\Uparrow }\right) \right) . \end{aligned}$$

  (57)

The time advancing formula is then limited to a stencil of only three cells ${\fancyscript{C}}_{i-1}$, ${\fancyscript{C}}_{i}$ and ${\fancyscript{C}}_{i+1}$, but the computation of the slopes for the cells ${\fancyscript{C}}_{i-1}$ and ${\fancyscript{C}}_{i+1}$ requires (see (53)) also to know the solution on the two sourrounding cells ${\fancyscript{C}}_{i-2}$ and ${\fancyscript{C}}_{i+2}$. The average solution $\bar{u}_{i}^{n+1}$, on each cell ${\fancyscript{C}}_i$ at time $t_{n+1} = t_n + \varDelta t$, can be computed knowing the solution on the augmented stencil, the $\mathrm {PV}$ introduced in Sect. 4, $\left\{ \bar{u}_{i-2}^n, \bar{u}_{i-1}^n, \bar{u}_{i}^n, \bar{u}_{i+1}^n, \bar{u}_{i+2}^n \right\}$ which constitutes the stencil obtained via the physical assembling algorithm (see Algorithm 3).