Dynamical low-rank approximation for Burgers' equation with uncertainty

05/10/2021 ∙ by Jonas Kusch, et al. ∙ KIT Leopold Franzens Universität Innsbruck Universität Tübingen 0

Quantifying uncertainties in hyperbolic equations is a source of several challenges. First, the solution forms shocks leading to oscillatory behaviour in the numerical approximation of the solution. Second, the number of unknowns required for an effective discretization of the solution grows exponentially with the dimension of the uncertainties, yielding high computational costs and large memory requirements. An efficient representation of the solution via adequate basis functions permits to tackle these difficulties. The generalized polynomial chaos (gPC) polynomials allow such an efficient representation when the distribution of the uncertainties is known. These distributions are usually only available for input uncertainties such as initial conditions, therefore the efficiency of this ansatz can get lost during runtime. In this paper, we make use of the dynamical low-rank approximation (DLRA) to obtain a memory-wise efficient solution approximation on a lower dimensional manifold. We investigate the use of the matrix projector-splitting integrator and the unconventional integrator for dynamical low-rank approximation, deriving separate time evolution equations for the spatial and uncertain basis functions, respectively. This guarantees an efficient approximation of the solution even if the underlying probability distributions change over time. Furthermore, filters to mitigate the appearance of spurious oscillations are implemented, and a strategy to enforce boundary conditions is introduced. The proposed methodology is analyzed for Burgers' equation equipped with uncertain initial values represented by a two-dimensional random vector. The numerical results show a reduction of the memory requirements, and that the important characteristics of the original system are well captured.



There are no comments yet.


page 19

page 20

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

A vast amount of engineering applications such as hydrology, gas dynamics, or radiative transport are governed by hyperbolic conservation laws. In many applications of interest the inputs (e.g. initial-, boundary conditions, or modeling parameters) of these equations are uncertain. These uncertainties arise from modeling assumptions as well as measurement or discretization errors, and they can heavily affect the behavior of inspected systems. Therefore, one core ingredient to obtain reliable knowledge of a given application is to derive methods which include the effects of uncertainties in the numerical simulations.

Methods to quantify effects of uncertainties can be divided into intrusive and non-intrusive methods. Non-intrusive methods run a given deterministic solver for different realizations of the input in a black-box manner, see e.g. xiu2005high; babuvska2007stochastic; loeven2008probabilistic for the stochastic-Collocation method and heinrich2001multilevel; mishra2012multi; mishra2012sparse; mishra2016numerical for (multi-level) Monte Carlo methods. Intrusive methods perform a modal discretization of the solution and derive time evolution equations for the corresponding expansion coefficients. The perhaps most prominent intrusive approach is the stochastic-Galerkin (SG) method ghanem2003stochastic, which represents the random dimension with the help of polynomials. These polynomials, which are picked according to a chosen probability distribution, are called generalized polynomial chaos (gPC) functions wiener1938homogeneous; xiu2002wiener

. By performing a Galerkin projection, a set of deterministic evolution equations for the expansion coefficients called the moment system can be derived.

Quantifying uncertainties in hyperbolic problems comes with a large number of challenges such as spurious oscillations le2004uncertainty; barth2013non; dwight2013adaptive or the loss of hyperbolicity poette2009uncertainty. A detailed discussion and numerical comparison of these challenges when using intrusive and non-intrusive methods can be found in kusch2020intrusive. Intrusive methods which preserve hyperbolicity are the intrusive polynomial moment (IPM) method poette2009uncertainty which performs a gPC expansion of the entropy variables, and the Roe transformation method pettersson2014stochastic; gerster2020entropies which performs a gPC expansion of the Roe variables. Furthermore, admissibility of the solution can be achieved by using bound-preserving limiters schlachter2018hyperbolicity to push the solution into an admissible set. Oscillations that frequently arise in hyperbolic problems can be mitigated by either filters kusch2020filtered or the use of multi-elements wan2006multi; durrwachter2020hyperbolicity; kusch2020oscillation.

One key challenge of uncertainty quantification is the exponential growth of the number of unknowns when the dimension of the random domain increases. Hyperbolic conservation laws, which tend to form shocks amplify this effect, since they require a fine discretization in each dimension. This does not only yield higher computational costs, but also extends memory requirements. Therefore, a crucial task is to find an efficient representation which yields a small error when using a small number of expansion coefficients. The gPC expansion provides such an efficient representation if the chosen polynomials belong to the probability density of the solution xiu2002wiener. However, the probability density is commonly known only for the initial condition and the time evolution of this density is not captured by the chosen gPC polynomials. Thus, choosing gPC polynomials according to the initial distribution can become inaccurate when the distribution changes over time.

We aim to resolve these issues by applying the dynamical low-rank approximation (DLRA) koch2007dynamical to our problem. To decrease memory and computational requirements DLRA represents and propagates the solution in time on a prescribed low-rank manifold. Such a low-rank representation is expected to be efficient, since choosing the gPC expansion according to the underlying probability density provides an accurate solution representation for a small number of expansion coefficients xiu2002wiener. Solving the DLRA equation by a matrix projector-splitting integrator lubich2014projector yields time evolution equations for updating spatial and uncertain basis functions in time. Hence, the resulting method is able to automatically adapt basis function to resolve important characteristics of the inspected problem. Furthermore, the matrix projector-splitting integrator has improved stability properties and error bounds kieri2016discretized. An extension of the matrix projector-splitting integrator to function spaces is presented in einkemmer2018low. First applications of dynamical low-rank approximations to uncertainty quantification are sapsis2009dynamically; sapsis2012dynamical; ueckermann2013numerical which use a so-called dynamical double orthogonal (DDO) approximation. The method however depends on the regularity of the coefficient matrix, which can potentially restrict the time step. Applications of the DLRA method in combination with the matrix projector-splitting integrator for parabolic equations with uncertainty can for example be found in musharbash2015error; musharbash2018dual; kazashi2020stability; kazashi2020existence. A dynamical low-rank approximation for random wave equations has been studied in musharbash2017symplectic. Examples of DLRA for kinetic equations, which include hyperbolic advection terms are einkemmer2018low; einkemmer2019quasi; einkemmer2020low; einkemmer2021asymptotic; peng2019low; peng2020high. Similar to our work, filters have been used in peng2019low to mitigate oscillations in the low-rank approximation. Furthermore, einkemmer2020low uses diffusion terms to dampen oscillatory artifacts.

In this work, we focus on efficiently applying the matrix projector-splitting integrator and the unconventional integrator for DLRA to hyperbolic problems with quadratic physical fluxes. Furthermore, we study the effects of oscillatory solution artifacts, which we aim to mitigate through filters. Furthermore, we investigate a strategy to preserve boundary conditions, similar to musharbash2018dual. Additionally, we investigate the unconventional DLRA integrator ceruti2020unconventional in the context of uncertainty quantification and compare it to the standard matrix projector-splitting integrator lubich2014projector

. The different integrators for dynamical low-rank are compared to the classical stochastic-Galerkin method for Burgers’ equation with two-dimensional uncertainties. In our numerical experiments, the chosen methods for DLRA capture the highly-resolved stochastic-Galerkin results nicely. It is observed that the unconventional integrator smears out the solution, which improves the expected value of the approximation but yields heavy dampening for the variance.

Following the introduction, we briefly present the required background for this paper in Section 2. Here, we give an overview of intrusive UQ methods as well as the DLRA framework applied to uncertainty quantification. Section 3 discusses the matrix projector-splitting integrator applied to a scalar, hyperbolic equation with uncertainty. Section 4 proposes a strategy to enforce Dirichlet boundary conditions in the low-rank numerical approximation and Section 5 discusses the numerical discretization. In Section 6 we demonstrate the effectiveness of the dynamical low-rank approximation ansatz for hyperbolic problems by investigating Burgers’ equation with a two-dimensional uncertainty.

2 Background

We compute a low-rank solution of a scalar hyperbolic equation with uncertainty


The solution depends on time , space

and a scalar random variable

. The random variable

is equipped with a known probability density function


2.1 Intrusive methods for uncertainty quantification

The core idea of most intrusive methods is to represent the solution to (1) by a truncated gPC expansion


Here, the basis functions are chosen to be orthonormal with respect to the probability density function , i.e.,

Note that the chosen polynomial ansatz yields an efficient evaluation of quantities of interest such as expected value and variance

This ansatz is used to represent the solution in (1) which yields


Then, a Galerkin projection is performed. We project the resulting residual to zero by multiplying (3) with test functions (for ) and taking the expected value. Prescribing the residual to be orthogonal to the test functions gives

This closed deterministic system is called the stochastic-Galerkin (SG) moment system. It can be solved with standard finite volume or discontinuous Galerkin methods, provided it is hyperbolic. The resulting moments can then be used to evaluate the expected solution and its variance. However, hyperbolicity is not guaranteed for non-scalar problems, which is why a generalization of SG has been proposed in poette2009uncertainty. This generalization, which is called the intrusive polynomial moment (IPM) method performs the gPC expansion on the so-called entropy variable , where is a convex entropy to (1). For more details on the IPM method, we refer to the original IPM paper poette2009uncertainty as well as (poette2019contribution, Chapter 4.2.3) and (10.5445/IR/1000121168, Chapter 1.4.5). Furthermore, the solutions of various UQ methods, including stochastic-Galerkin, show spurious solution artifacts such as non-physical step wise approximations le2004uncertainty; barth2013non; dwight2013adaptive. One strategty to mitigate these effects are filters, which have been proposed in kusch2020filtered. The idea of filters is to dampen high order moments in between time steps to ensure a smooth approximation in the uncertain domain .

When denoting the moments evaluated in a spatial cell at time as , a finite volume update with numerical flux takes the form


To reduce computational costs, one commonly chooses a kinetic flux

where is a numerical flux for the original problem (1) and we make use of a quadrature rule with points and weights . An efficient computation is achieved by precomputing and storing the terms . For a more detailed derivation of the finite volume scheme, see e.g. kusch2020filtered. Now, to mitigate spurious oscillations in the uncertain domain, we derive a filtering matrix . Note that to approximate a given function , the common gPC expansion minimizes the L-distance between the polynomial approximation and . Unfortunately, this approximation tends to oscillate. A representation

which dampens oscillations in the polynomial approximation can be derived by solving the optimization problem


Here, is a user-determined filter strength and is an operator which returns high values when the approximation oscillates. Note that we added a term which punishes oscillations in the L

-distance minimization that is used in gPC expansions. For uniform distributions, a common choice of


In this case, the gPC polynomials are eigenfunctions of

. In this case, the optimal expansion coefficients are given by , i.e., high order moments of the function will be dampened. Collecting the dampening factors in a matrix and applying this dampening steps in between finite volume updates yields the filtered SG method


A filter for high dimensional uncertainties can be applied by successively adding a punishing term to the optimization problem (5) for every random dimension. When the random domain is two-dimensional, i.e., , then with

one can use

Here, the gPC functions are tensorized and collected in a vector

where . The resulting filtered expansion coefficients are then given by

The costs of the filtered SG as well as the classical SG method when using a kinetic flux function are and the memory requirement is . For more general problems with uncertain dimension , tensorizing the chosen basis functions and quadrature yields and the memory requirement is . For a discussion of the advantages when using a kinetic flux, see (kusch2020intrusive, Appendix A).

2.2 Matrix projector-splitting integrator for dynamical low-rank approximation

The core idea of the dynamical low-rank approach is to project the original problem on a prescribed manifold of rank functions. Such an approximation is given by


In the following, we denote the set of functions, which have a representation of the form (7) by . Then, instead of computing the best approximation in the aim is to find a solution fulfilling


Here, denotes the tangent space of at . According to einkemmer2018low, the orthogonal projection onto the tangent space reads


This lets us rewrite (8) as


where denotes the integration over the spatial domain. A detailed derivation can be found in (koch2007dynamical, Lemma 4.1). Then, a Lie-Trotter splitting technique yields


In these split equations, each solution has a decomposition of the form (7). The solution of these split equations can be further simplified: Let us write , i.e., we define . Then, we test (10a) against and omit the index in the decomposition to simplify notation. This gives


This system is similar to the SG moment system, but with time dependent basis functions, cf. tryoen2010intrusive. Performing a Gram-Schmidt decomposition in of yields time updated and . Then testing (10b) with and gives


This equation is used to update . Lastly, we write and test (10c) with . This then yields


Again, a Gram-Schmidt decomposition is used on to determine the time updated quantities and . Note that the -step does not modify the basis and the -step does not modify the basis . Furthermore, the -step solely alters the coefficient matrix . Therefore, the derived equations can be interpreted as time update equations for the spatial basis functions , the uncertain basis functions and the expansion coefficients . Hence, the matrix projector-splitting integrator evolves the basis functions in time, such that a low-rank solution representation is maintained. Let us use Einstein’s sum notation to obtain a compact presentation of the integrator. The projector splitting procedure which updates the basis functions , and coefficients from time to then takes the following form:

  1. -step: Update to and to via

    Determine and with .

  2. -step: Update via

    and set .

  3. -step: Update and via

    Determine and with .

2.3 Unconventional integrator for dynamical low-rank approximation

Note that the -step (12) evolves the matrix backward in time, which is a source of instability for non-reversible problems such as diffusion equations or particle transport with high scattering. Furthermore, the presented equation must be solved successively, which removes the possibility of solving different steps in parallel. In ceruti2020unconventional, a new integrator which enables parallel treatment of the and -step while only evolving the solution forward in time has been proposed. This integrator, which is called the unconventional integrator, is similar but not equal to the matrix projector-splitting integrator and works as follows:

  1. -step: Update to via


    with a QR-decomposition

    and store .

  2. -step: Update to via

    Determine with a QR-decomposition and store .

  3. -step: Update to via

    and set .

Note that the first and second steps can be performed in parallel. Furthermore, the unconventional integrator inherits the exactness and robustness properties of the classical matrix projector-splitting integrator, see ceruti2020unconventional. Additionally, it allows for an efficient use of rank adaptivity CKL2021.

3 Matrix projector-splitting integrator for uncertain hyperbolic problems

Before deriving the evolution equations for scalar problems, we first discuss hyperbolicity.

Theorem 1.

Provided the original scalar equation (1) is hyperbolic, its corresponding -step equation (11

) is hyperbolic as well, i.e., its flux Jacobian is diagonalizable with real eigenvalues.


First, we apply the chain rule to the

-step equation to obtain

Linearity of the expected value gives

Then, the flux Jacobian is obviously symmetric, i.e., by the spectral theorem it is diagonalizable with real eigenvalues. ∎

We remind the reader that the -step of the unconventional integrator equals the -step of the matrix projector-splitting integrator. Therefore, Theorem 1 holds for both numerical integrators. Note that this result only holds for scalar equations. In the system case, hyperbolicity cannot be guaranteed. Methods to efficiently guarantee hyperbolicity for systems are not within the scope of this paper and will be left for future work.

Now, we apply the matrix projector-splitting integrator presented in Section 2.2 to Burgers’ equation


Our goal is to never compute the full solutions , and but to rather work on the decomposed quantities saving memory and reducing computational costs. Here, we represent the uncertain basis functions with gPC polynomials, i.e.,


The spatial basis functions are represented with the help of finite volume (or DG0) basis functions


where is the indicator function on the interval . Then, the spatial basis functions are given by


In the following, we present an efficient evaluation of non-linear terms that arise in the matrix projector-splitting integrator for the non-linear equations.

3.1 -step

Let us start with the -step, which decomposes the solution into

Plugging this representation as well as the quadratic flux of Burgers’ equation into (11) gives

Defining yields the differential equation


which according to Theorem 1 is guaranteed to be hyperbolic. Representing the spatial coordinate with points, a number of evaluations per time step is needed for the evaluation of (18). The terms can be computed by applying a quadrature rule with weights and points . Then, we can compute for and in operations. Furthermore, one can compute


in operations. Hence, the total costs for the -step, which we denote by are

It is important to point out the essentially quadratic costs of . This non-linear term stems from the modal gPC approximation of the uncertain basis . We will discuss a nodal approximation in Section 3.5, which yields linear costs. If we denote the memory requires for the -step by , we have

when precomputing and storing the terms .

3.2 -step

For the -step, the basis functions and remain constant in time, i.e., we have


Plugging this representation into the -step (12) yields


Using a sufficiently accurate quadrature rule enables an efficient computation of the terms . The values of for and can be computed in operations, since

The terms can be precomputed and stored. Furthermore, we make use of . Then, a sufficiently accurate quadrature (e.g. a Gauss quadrature with ) yields

in operations. The derivative term in can be approximated with a finite volume stencil or a simple finite difference stencil. Dividing the spatial domain into points lets us define elements . Then, at spatial position we choose the finite difference approximation


Again choosing DG0 basis functions yields

Summing over spatial cells to approximate the spatial integral gives


The required number of operations to compute this term for is . Furthermore, the memory requirement is to store all as well as to store all integral terms (23). Then the multiplication in (21) requires operations. Hence the total costs and memory requirements for the -step are

3.3 -step

The final step is the -step (13), which for Burgers’ equation becomes


To obtain a time evolution equation of the expansion coefficients such that


we test (3.3) with , which gives


The terms can be computed analogously to (19) in operations, taking up memory of . Precomputing the terms again has memory requirements of . The term can be reused from the -step computation and the multiplication in (26) requires operations. Therefore, the overall costs and memory requirements for the -step are

Note that the QR decompositions needed to compute from and require as well as operations respectively in every time step.

3.4 Filtered Matrix projector-splitting integrator

Similar to filtered stochastic-Galerkin, we wish to apply a filtering step in between time steps. Let us write the low-rank solution as

Following the idea of filtering, we now wish to determine such that the solution representation minimizes the L-error in combination with a term punishing oscillatory solution values. Equivalently to the derivation of the fSG ansatz (5), this gives

As before we have