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.
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 ofis
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 vectorwhere . 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:
-step: Update to and to via
Determine and with .
-step: Update via
and set .
-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:
-step: Update to via
with a QR-decompositionand store .
-step: Update to via
Determine with a QR-decomposition and store .
-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.
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.
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 .
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
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