Oscillation Mitigation of Hyperbolicity-Preserving Intrusive Uncertainty Quantification Methods for Systems of Conservation Laws
OOscillation Mitigation of Hyperbolicity-Preserving Intrusive UncertaintyQuantification Methods for Systems of Conservation Laws
Jonas Kusch a , Louisa Schlachter b a Computational Science and Mathematical Methods, Karlsruhe Institute of Technology, Englerstraße 2, 76128 Karlsruhe, [email protected] b Fachbereich Mathematik, TU Kaiserslautern, Erwin-Schr¨odinger-Str., 67663 Kaiserslautern, Germany, [email protected]
Abstract
In this article we study intrusive uncertainty quantification schemes for systems of conservation laws withuncertainty. Standard intrusive methods lead to oscillatory solutions which sometimes even cause the lossof hyperbolicity. We consider the stochastic Galerkin scheme, in which we filter the coefficients of thepolynomial expansion in order to reduce oscillations. We further apply the multi-element approach andensure the preservation of hyperbolic solutions through the hyperbolicity limiter. In addition to that,we study the intrusive polynomial moment method, which guarantees hyperbolicity at the cost of solvingan optimization problem in every spatial cell and every time step. To reduce numerical costs, we applythe multi-element ansatz to IPM. This ansatz decouples the optimization problems of all multi elements.Thus, we are able to significantly decrease computational costs while improving parallelizability. We finallyevaluate these oscillation mitigating approaches on various numerical examples such as a NACA airfoil anda nozzle test case for the two-dimensional Euler equations. In our numerical experiments, we observe themitigation of spurious artifacts. Furthermore, using the multi-element ansatz for IPM significantly reducescomputational costs.
Keywords:
Uncertainty quantification, intrusive, polynomial Chaos, stochastic Galerkin, Filter, IPM,Hyperbolicity, Limiter, Multi-Element
1. Introduction
Uncertainties play a key role in modelling hyperbolic systems of conservation laws if crucial data or param-eters might not be available exactly due to measurement errors. Thus, to achieve an adequate descriptionof reality, one needs to include non-deterministic effects in the approximation of deterministic systems. Inthe recent years many approaches, e.g., methods based on Bayesian inversion, Monte Carlo algorithms orstochastic Galerkin schemes, were proposed to quantify the uncertainties in order to account for them inpredictions and simulations [1, 2, 15, 22, 23, 49, 55, 64]. In the context of these Uncertainty Quantification(UQ) methods, we distinguish between so-called non-intrusive and intrusive approaches.Non-intrusive UQ methods use deterministic solvers of the equations as a black box code in order to solve themodel problem by repeated application of this black box on fixed realizations of the uncertainty. The mostwidely known non-intrusive UQ method is (Multi-Level) Monte Carlo [14, 19, 33] which is based on statisticalsampling methods and can be easily implemented and adopted to any type of uncertain conservation law,but comes with potentially high cost due to repeated application of e.g. finite volume methods (FVM).These so-called MC-FVM schemes have been studied in [37, 40] showing a slow convergence rate that is
Preprint submitted to
Journal of Computational Physics a r X i v : . [ m a t h . NA ] A ug mproved by Multi-Level MC-FVM algorithms for conservation laws [38, 39]. Another non-intrusive UQscheme is stochastic collocation [62], further developments of this method are described in [34, 41, 61].In this article, we focus on intrusive UQ methods, which aim to increase the overall efficiency but requirea modification of the underlying solver for the deterministic problem. The most popular methods are thestochastic Galerkin (SG) and the intrusive polynomial moment (IPM) method [46]. They both rely on thegeneralized Polynomial Chaos (gPC) expansion [1, 9, 57, 63] which is theoretically based on the PolynomialChaos expansion from [60]. Stochastic Galerkin expands the solution in the conserved variables such thatthe gPC system results in a weak formulation of the equations with respect to the stochastic variable. IPMexpands the stochastic solution in the so-called entropic variables, which results in a hyperbolic gPC systemthat yields a good approximation quality, however, it is necessary to know a strictly convex entropy of theoriginal system beforehand.The biggest challenge of UQ methods for hyperbolic conservation laws lies in the fact that discontinuities inthe physical space propagate into the solution manifold [6] and cause oscillations. The naive usage of SG fornonlinear hyperbolic problems even typically fails [3, 46], since the polynomial expansion of discontinuousdata yields huge Gibbs oscillations that result in the loss of hyperbolicity. In order to resolve this problem, weapply the so-called hyperbolicity-preserving limiter to the classical SG approach. This limiter is introducedin [50] and the theory is extended in [13] to high-order Runge-Kutta discontinuous Galerkin schemes aswell as the so-called multi-element approach developed by [53, 58], where the random space is divided intodisjoint elements in order to define local gPC approximations. The multi-element ansatz corresponds toan h -refinement in the random space. Further developments of this method encompass h - and hp -adaptiverefinements in the stochastic space [53, 56, 59]. As presented in [51], a simple example of the uncertainlinear transport equation shows that the multi-element ansatz and the hyperbolicity limiter is still notenough to sufficiently reduce oscillations. Similar situations have been observed in [5, 34, 46]. Anotherapproach to damp oscillations induced by the Gibbs phenomenon is given in [29], where filters are appliedto the gPC coefficients of the SG approximation. Filters are a common technique from kinetic theory [35]to reduce oscillations but do not guarantee the hyperbolicity of the gPC system. We therefore employ themulti-element ansatz together with the hyperbolicity limiter and the filter to stochastic Galerkin in orderto obtain a robust numerical scheme that is able to deal with Gibbs phenomenon.On the other hand, the hyperbolicity of the moment system is ensured through the IPM method that has beenintroduced in [46]. Expanding the stochastic solution not in the conserved variables but in entropic variablesis well known in the radiative transfer community as minimum entropy models [17, 27, 32]. The resulting gPCsystem has good approximation properties but requires to solve typically expensive optimization problemsfor nonlinear systems in every space-time cell in order to calculate the entropic variable. IPM may beused to control oscillations of the solution since they are bounded to a certain range though the entropy ofthe system. Another development of this method is proposed in [26], where a second-order IPM scheme isdescribed which fulfills the maximum principle. More information on IPM can be found in [12, 27, 28, 31, 47].Within the scope of this work, we adapt the multi-element ansatz on the IPM scheme. As described above,this approach is already successfully used within the SG system. Dividing the stochastic domain into multi-elements promises to reduce oscillations and to decrease computational times within IPM since it decouplesand simplifies the optimization problems. In order to further reduce the numerical costs we additionallyapply the one-shot implementation for IPM from [31], similar to the idea of one-shot optimization in shapeoptimization [18]. Moreover, each of our numerical schemes is accelerated employing adaptivity, i.e. thetruncation order of the gPC polynomials adapts locally to the smoothness of the solution, cf. [24, 31, 36, 54].The paper is structured as follows. In Section 2 we describe our problem setting for systems of hyperbolicconservation laws with uncertainty. We then describe in Section 3 the filtered hyperbolicity-preservingstochastic Galerkin scheme by introducing the multi-element approach as well as the hyperbolicity limiterand filters for stochastic Galerkin. In Section 4 we complete the theoretical framework by deriving the multi-element intrusive polynomial moment method. With these oscillation mitigating intrusive UQ schemes athand, we present numerical test cases for the one- and two-dimensional Euler equations within Section 5.2ere, we study the behavior of the proposed methods in terms of errors and computational costs, thatverify the reduction of oscillations compared to the standard intrusive schemes. Moreover, we are able tosignificantly reduce the numerical costs by employing the multi-element ansatz in IPM.
2. Problem Setting
We consider random systems of hyperbolic conservation laws of the form ∂ t u ( t, x , ξ ) + ∇ x · F ( u ( t, x , ξ ) , ξ ) = , (2.1a)for a flux function F = ( f , . . . , f M ), with f m ( u ) : R d → R d , m = 1 , . . . , M , for x ∈ X ⊂ R M , and where thesolution u = u ( t, x , ξ ) : R + × R M × Ξ → R d is depending on a one-dimensional random variable ξ with probability space (Ω , F , P ). Here, we denote therandom space of this uncertainty by Ξ := ξ (Ω) and its probability density function by f Ξ ( ξ ) : Ξ → R + . Wefurther assume that the uncertainty can be introduced via the initial conditions, namely u ( t = 0 , x , ξ ) = u ( x , ξ ) . (2.1b)Depending on X , additional boundary conditions have to prescribed. In this article, we consider intrusiveuncertainty quantification methods in order to solve (2.1). Therefore, we describe the most commonly usedstochastic Galerkin and intrusive polynomial moment method in the following sections, where we modifythem according to the theory explained within the introduction. The proposed methods are acceleratedusing techniques discussed in [31].
3. Filtered hyperbolicity-preserving stochastic Galerkin scheme
In this section we explain oscillation mitigating approaches for the stochastic Galerkin scheme. As mentionedin Section 1, stochastic Galerkin systems are likely to lose hyperbolicity [46], since discontinuities thattypically arise in hyperbolic equations propagate into the stochastic domain and cause Gibbs oscillationsthat might yield the loss of hyperbolicity, i.e. the failure of the standard stochastic Galerkin scheme. Wetherefore apply the multi-element approach [53, 58] as well as the hyperbolicity-preserving limiter thatwas introduced in [50]. However, this limiter still does not sufficiently damp Gibbs oscillations as seenin [51], where an uncertain linear advection example reveals huge overshoots that are not reduced by thehyperbolicity limiter due to the preset hyperbolicity in this scalar equation. In this article, we consideranother approach that aims on reducing oscillations, i.e. we combine the filtered SG method [29] with thehyperbolicity limiter such that the method can be applied on any system of conservation laws. To thisend, we filter the gPC coefficients of the SG approximation, which allows a numerically cheap reduction ofoscillations. Different filters for SG schemes are derived and described in [29].
We seek for an approximate solution by a finite-term generalized Polynomial Chaos expansion (see e.g. [15]) u ( t, x , ξ ) ≈ K Ξ (cid:88) k =0 u k ( t, x ) φ k ( ξ ) , (3.1)where the polynomials φ k of degree k are supposed to satisfy the orthogonality relation (cid:90) Ξ φ k ( ξ ) φ ˜ k ( ξ ) f Ξ ( ξ )d ξ = (cid:40) , if k = ˜ k , else ∀ k, ˜ k ∈ { , . . . , K Ξ } . (3.2)3nserting (3.1) into (2.1) and applying a Galerkin projection in the stochastic space leads to the so-calledstochastic Galerkin system ∂ t u ˜ k ( t, x ) + ∇ x · (cid:90) Ξ F (cid:32) K Ξ (cid:88) k =0 u k ( t, x ) φ k ( ξ ) , ξ (cid:33) φ ˜ k ( ξ ) f Ξ ( ξ )d ξ = 0 , ˜ k = , . . . , K Ξ . (3.3)For discontinuous solutions, the gPC approach may converge slowly or even fail to converge, cf. [46, 56].As presented in [13, 58], we therefore apply the multi-element approach, where Ξ is divided into disjointelements with local gPC approximations of (2.1).We assume that Ξ = ( ξ L , ξ R ) and define a decomposition of Ξ into N Ξ multi-elements Ξ l = ( ξ l − , ξ l + ), l = 1 , . . . , N Ξ , of width ∆ ξ = ξ L − ξ R N Ξ . Moreover, we introduce an indicator variable χ l : Ω → { , } on everyrandom element χ l ( ξ ) := (cid:40) ξ ∈ Ξ l , l = 1 , . . . , N Ξ . If we let { φ k,l ( ξ ) } ∞ k =0 be orthonormal polynomials with respect to a conditional probabilitydensity function f Ξ l : Ξ l → R on the multi-element Ξ l , as in [13], the global approximation (3.1) can bewritten as u ( t, x , ξ ) = N Ξ (cid:88) l =1 u l ( t, x , ξ ) χ l ( ξ ) ≈ N Ξ (cid:88) l =1 K Ξ (cid:88) k =0 u k,l ( t, x ) φ k,l ( ξ ) χ l ( ξ ) . (3.5)As N Ξ , K Ξ → ∞ , the local approximation converges to the global solution in L (Ω), cf. [4]. Remark 3.1.
For a uniform distribution of the uncertainty, we have f Ξ l ( ξ ) = ξ l +1 − ξ l χ Ξ l ( ξ ) . Moreover,the local probability density function remains a density function of a uniformly distributed random variable.Other distribution types such as the normal distribution require to numerically compute a set of polynomialswhich are orthogonal with respect to f Ξ l ( ξ ) , see [58]. We come up with the following ME stochastic Galerkin system ∂ t u ˜ k,l ( t, x ) + ∇ x · (cid:90) Ξ l F (cid:16) N Ξ (cid:88) l =1 K Ξ (cid:88) k =0 u k,l ( t, x ) φ k,l ( ξ ) , ξ (cid:17) φ ˜ k,l ( ξ ) f Ξ l ( ξ )d ξ = , (3.6)for ˜ k = 0 , . . . , K Ξ and l = 1 , . . . , N Ξ .The calculation of the expected value and variance of u can be found in [13]. We now describe our numerical scheme, whereas we consider two-dimensional physical spaces which werequire for our numerical examples of the two-dimensional Euler equations. Thus, we divide the spatialdomain X = [ x L , x R ] × [ y L , y R ] ⊂ R into a uniform rectangular mesh with cells X i,j = [ x i − , x i + ] × [ y j − , y j + ] where x i ± = x i ± ∆ x , y j ± = y j ± ∆ y and ∆ x = x R − x L N x , ∆ y = y R − x L N y , such that ( i, j ) ∈ N X := { ( i, j ) ∈ N | i ≤ N x , j ≤ N y } .Additionally, we discretize the random space Ξ with the multi-element ansatz from Section 3.1, i.e. wedivide Ξ into Ξ l , for l = 1 , . . . , N Ξ . We test (2.1) in each spatial cell X i,j by a test function v ( x, y ) with4upp( v ) ⊆ X i,j and obtain after one formal integration by parts the following weak formulation ∂ t (cid:90) X i,j u ( t, x, ξ ) v ( x, y )d( x, y ) − (cid:90) X i,j (cid:16) f (cid:0) u ( t, x, ξ ) , ξ (cid:1) ∂ x v ( x, y ) + f (cid:0) u ( t, x, ξ ) , ξ (cid:1) ∂ y v ( x, y ) (cid:17) d( x, y ) (3.7a)+ y j + 12 (cid:90) y j − (cid:104) f (cid:0) u ( t, x, ξ ) , ξ (cid:1) v ( x, y ) (cid:105) x i + 12 x i − d y (3.7b)+ x j + 12 (cid:90) x j − (cid:104) f (cid:0) u ( t, x, ξ ) , ξ (cid:1) v ( x, y ) (cid:105) y i + 12 y i − d x = , (3.7c)for ( i, j ) ∈ N X .Since the solution is discontinuous, we replace the flux at the interfaces x i ± with a consistent numericalflux function ˆ f m , such that ˆ f m ( u , u , ξ ) = f m ( u , ξ ) for m = 1 ,
2. If we use a first order numerical scheme,we choose v = x ∆ y and use the midpoint rule to approximate the integrals in (3.7b) and (3.7c). We thenobtain for one time step of the forward Euler scheme u ( n +1) i,j = u ( n ) i,j − ∆ t ∆ x (cid:16) ˆ f ( u ( n ) i,j , u ( n ) i +1 ,j , ξ ) − ˆ f ( u ( n ) i − ,j , u ( n ) i,j , ξ ) (cid:17) (3.8a) − ∆ t ∆ y (cid:16) ˆ f ( u ( n ) i,j , u ( n ) i,j +1 , ξ ) − ˆ f ( u ( n ) i,j − , u ( n ) i,j , ξ ) (cid:17) , (3.8b)where u ( n ) i,j , ( i, j ) ∈ N X denotes the spatial cell mean of u in X i,j at time t n , namely u ( n ) i,j = u ( n ) i,j ( ξ ) := 1∆ x ∆ y (cid:90) X i,j u ( t n , x, y, ξ ) d( x, y ) . In our numerical experiments we use the numerical Lax-Friedrichs fluxˆ f ( u i,j , u i +1 ,j , ξ ) = 12 (cid:0) f ( u i,j , ξ ) + f ( u i +1 ,j , ξ ) − λ max1 ( u i +1 ,j − u i,j ) (cid:1) , ˆ f ( u i,j , u i,j +1 , ξ ) = 12 (cid:0) f ( u i,j , ξ ) + f ( u i,j +1 , ξ ) − λ max2 ( u i,j +1 − u i,j ) (cid:1) , where the numerical viscosity constants λ max1 and λ max2 are taken as the global estimates of the absolutevalue of the largest eigenvalue of ∂ f ( u ) ∂ u and ∂ f ( u ) ∂ u , respectively. The CFL condition of (3.8) reads∆ t ≤ ∆ xλ max1 + ∆ yλ max2 . For the coefficients of the gPC polynomial (3.5) in the multi-element stochastic Galerkin system (3.6), thefinite volume scheme (3.8) is given by u ( n +1) k,i,j,l = u ( n ) k,i,j,l − ∆ t ∆ x (cid:90) Ξ l (cid:16) ˆ f ( u ( n ) i,j , u ( n ) i +1 ,j , ξ ) − ˆ f ( u ( n ) i − ,j , u ( n ) i,j , ξ ) (cid:17) φ k,l f Ξ l d ξ (3.9a) − ∆ t ∆ y (cid:90) Ξ l (cid:16) ˆ f ( u ( n ) i,j , u ( n ) i,j +1 , ξ ) − ˆ f ( u ( n ) i,j − , u ( n ) i,j , ξ ) (cid:17) φ k,l f Ξ l d ξ, (3.9b)5or each equation k = 0 , . . . , K Ξ , where u ( n ) k,i,j,l describes the X i,j × Ξ l cell mean of the k th coefficient in thegPC expansion at time t n , namely u ( n ) k,i,j,l := 1∆ x ∆ y (cid:90) X i,j (cid:90) Ξ l u ( t n , x, ξ ) φ k,l ( ξ ) f Ξ l ( ξ )d ξ d( x, y ) , for k = 0 , . . . , K Ξ , ( i, j ) ∈ N X and l = 1 , . . . , N Ξ . Remark 3.2.
The integrals in (3.9) are numerically solved using a Q Ξ -point appropriate quadrature ruleon Ξ . In our numerical examples we use Gauß-Legendre or Clenshaw-Curtis quadrature rules. Note thatcomputing moments of an underlying flux is a common numerical flux choice in kinetic theory [11, 16, 42, 43]and is frequently used in the context of uncertainty quantification, see e.g. [25, 50]. A comparison tonumerical fluxes working directly on the moments can be found in [10] or [31, Appendix A].3.3. Hyperbolicity-Preserving Limiter Usually, the solution of our system of equations (2.1) has to fulfill certain physical properties. For example, adensity should always be nonnegative. For hierarchical models like the Euler equations, the system normallyloses hyperbolicity for nonphysical states, cf. [46].
Definition 3.3.
We call the set R := (cid:40) u ∈ R d (cid:12)(cid:12) M (cid:88) m =1 α m ∂ f m ( u , ξ ) ∂ u has d real eigenvalues and is diagonalizable ∀ α , . . . , α M ∈ R (cid:41) the hyperbolicity set . We call every solution vector u ∈ R admissible . Assumption 3.4.
In the following we always assume that the hyperbolicity set R is open and convex. In order to apply the hyperbolicity-limiter from [50], we require positivity-preserving numerical fluxes ˆ f , ˆ f . Definition 3.5.
A scheme (3.8) and the numerical fluxes ˆ f , ˆ f are called positivity-preserving , if u ( n ) i,j ∈R for all ( i, j ) ∈ N X at time t n implies that u ( n +1) i,j ∈ R for all ( i, j ) ∈ N X at time t n +1 . Assumption 3.6.
The numerical fluxes ˆ f , ˆ f , are positivity-preserving under a suitable CFL condition λ max1 ∆ t ∆ x + λ max2 ∆ t ∆ y ≤ C, where C ∈ (0 , . Remark 3.7.
The positivity-preserving property of the Lax-Friedrichs numerical Flux has been proven in[52, 65]. For the Lax-Friedrichs numerical flux, one may choose C = 0 . , cf. [65, Remark 2.4] and [44]. According to [50], the numerical scheme (3.9) preserves hyperbolicity of the zeroth SG moment u ( n )0 ,i,j,l for positivity-preserving numerical fluxes and we can now define the hyperbolicity limiter which moves apossibly nonadmissible solution towards this admissible moment and inside the hyperbolicity set. Followingthe outline of [50], at time t n we define the slope-limited SG polynomial in the cell X i,j × Ξ l asΛΠ θK Ξ ( u ( n ) i,j ) (cid:12)(cid:12) Ξ l ( ξ ) := θ u ( n )0 ,i,j,l + (1 − θ ) K Ξ (cid:88) k =0 u ( n ) k,i,j,l φ k,l ( ξ )= u ( n )0 ,i,j,l + (1 − θ ) K Ξ (cid:88) k =1 u ( n ) k,i,j,l φ k,l ( ξ ) , (3.10)6or ( i, j ) ∈ N X and l = 1 , . . . , N Ξ . For more information see [50]. We chooseˆ θ ( n ) i,j,l := inf (cid:110) ˜ θ ∈ [0 , (cid:12)(cid:12) ΛΠ ˜ θK Ξ ( u ( n ) i,j ) (cid:12)(cid:12) Ξ l ∈ R for all ξ quadrature nodes (cid:111) , which ensures hyperbolicity of the slope-limited SG polynomial (3.10). Due to the openness of R , we needto modify θ slightly in order to avoid placing the solution onto the boundary (if the limiter was active).Therefore we use θ = (cid:40) ˆ θ, if ˆ θ = 0 , min(ˆ θ + (cid:15), , if ˆ θ > , where 0 < (cid:15) = 10 − should be chosen small enough to ensure that the approximation quality is notinfluenced significantly. Remark 3.8.
Calculating the slope-limited SG polynomial (3.10) , we can simply replace the SG momentsby ΛΠ θK Ξ ( u ( n ) k,i,j,l ) = (cid:40) u ( n )0 ,i,j,l , if k = 0 , (1 − θ ) u ( n ) k,i,j,l , if k > , k = 0 , . . . , K Ξ , where θ is chosen separately for each t n and X i,j × Ξ l , ensuring that ΛΠ θK Ξ ( u ( n ) i,j ) (cid:12)(cid:12) Ξ l ∈ R in all cells. Remark 3.9.
We may calculate the the value of the limiter variable θ directly by solving the minimizationproblem min θ (3.11a) s.t. θ ∈ [0 ,
1] (3.11b) θ (cid:101) u + (1 − θ ) u ∈ R , (3.11c) where (cid:101) u represents the admissible zeroth SG moment. An exemplary derivation of this parameter for thetwo-dimensional compressible Euler equations can be found in [13]. Example 3.10.
We consider the one-dimensional compressible Euler equations ∂ t ρρvρe + ∂ x ρvρv + pv ( ρe + p ) = . (3.12) These equations describe the evolution of a gas with density ρ , velocity v and specific total energy e . Thepressure of the gas can be determined from p = ( γ − ρ (cid:18) e − v (cid:19) and γ is the heat capacity ratio. The hyperbolicity set from Definition 3.3 is given by solutions with positivedensity and pressure, i.e. R = u = ρρvρe (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ρ > , p = ( γ − ρ (cid:18) e − v (cid:19) > . hen, the solution of the minimization problem (3.11) reads θ (cid:63) = max ( h ( θ ) , h ( θ ) , h ( θ − )) ,θ = ρρ − (cid:101) ρ ,θ ± = E (cid:101) ρ − E ρ + (cid:101) E ρ − m (cid:101) m + ( m ) (cid:0) m − (cid:101) m (cid:1) + 2 (cid:0) E (cid:101) ρ + (cid:101) E ρ − E ρ − (cid:101) E (cid:101) ρ (cid:1) ± (cid:113)(cid:0) E (cid:101) ρ − (cid:101) Eρ (cid:1) + 2 (cid:0) E ( (cid:101) m ) ρ + (cid:101) E ( m ) (cid:101) ρ − E m (cid:101) m (cid:101) ρ − (cid:101) E m (cid:101) m ρ (cid:1)(cid:0) m − (cid:101) m (cid:1) + 2 (cid:0) E (cid:101) ρ + (cid:101) E ρ − E ρ − (cid:101) E (cid:101) ρ (cid:1) ,h ( x ) = x [0 , ( x ) . The calculation of the above quantities follows similarly to [13].3.4. Filtered stochastic Galerkin
In addition to the hyperbolicity-preserving limiter from Section 3.3, we apply a filter to the stochasticGalerkin scheme, as introduced in [29]. Filters are a common technique in spectral methods [8, 20] andhave for example been applied in the context of kinetic theory [35] in order to dampen oscillations. Wewill similarly apply them to the gPC coefficients in (3.1) to reduce oscillations in the uncertainty, whilemaintaining high-order accuracy in deterministic regions.In the following we present the so-called L filter, adding a penalizing term to the L error of the solution inorder to reduce oscillations within high-order coefficients of the gPC expansion. We follow the outline from[29] and include it into the framework of the multi-element approach. The standard SG method chooses thecoefficients in the gPC expansion such that the cost functional J := 12 N Ξ (cid:88) l =1 (cid:90) Ξ l (cid:13)(cid:13) u l − u K Ξ l (cid:13)(cid:13) f Ξ l d ξ is minimal with respect to the multi-element gPC polynomial u K Ξ of degree K Ξ , where (cid:107) · (cid:107) denotes the L norm and u = u ( t, x, y, ξ ) with ξ ∈ Ξ l is the solution to the uncertain conservation law in the l thmulti-element for l = 1 , . . . , N Ξ . The minimizer is given by (3.5).The filtered stochastic Galerkin polynomial is now set to (cid:98) u K Ξ ( t, x, y, ξ ) := N Ξ (cid:88) l =1 K Ξ (cid:88) k =0 g ( k ) u k,l ( t, x, y ) φ k,l ( ξ ) χ l ( ξ ) , (3.13)with the filter function g . For the L filter, a filter function and therefore a filtered stochastic Galerkinsolution, is obtained by minimizing J λ := 12 N Ξ (cid:88) l =1 (cid:90) Ξ l (cid:13)(cid:13) u l − (cid:98) u K Ξ l (cid:13)(cid:13) f Ξ l d ξ + λ N Ξ (cid:88) l =1 (cid:90) Ξ l (cid:13)(cid:13) L ( (cid:98) u K Ξ l ) (cid:13)(cid:13) f Ξ l d ξ, (3.14)for the filter strength λ ≥
0. The operator L punishes oscillations and is chosen such that the basis poly-nomials of the gPC expansion are eigenfunctions of L . Differentiating (3.14) with respect to the coefficients (cid:98) u K Ξ k,l and setting the result equal to zero yields the filter function g .8 emark 3.11. For a uniform distribution of ξ , a common choice of the operator L is L ( u ( ξ )) = ((1 − ξ ) u (cid:48) ( ξ )) (cid:48) , since the Legendre polynomials are eigenfunctions of L , i.e. they fulfill L ( φ k,l ( ξ )) = − k ( k + 1) φ k,l ( ξ ) , for k = 0 , . . . , K Ξ and l = 1 , . . . , N Ξ . Other choices of the filter function are Lanczos and ErfcLog [48]. Wedifferentiate the cost functional (3.14) with respect to (cid:98) u K Ξ k,l and deduce d (cid:98) u K Ξ k,l J λ = (cid:90) Ξ l (cid:12)(cid:12)(cid:12)(cid:12) K Ξ (cid:88) ˜ k =0 ( u ˜ k,l − (cid:98) u K Ξ ˜ k,l ) φ ˜ k,l (cid:12)(cid:12)(cid:12)(cid:12) φ k,l f Ξ l d ξ − (cid:90) Ξ l (cid:12)(cid:12)(cid:12)(cid:12) K Ξ (cid:88) ˜ k =0 λk ( k + 1) (cid:98) u K Ξ ˜ k,l φ ˜ k,l (cid:12)(cid:12)(cid:12)(cid:12) φ k,l f Ξ l d ξ ! = . Hence, the filter function within the filtered SG polynomial (3.13) for a uniformly distributed uncertaintyreads g ( k ) = 11 + λk ( k + 1) . (3.15) Thus, the zeroth SG moment stays unfiltered due to g (0) = 1 which means that the cell means are preserved.For larger values of k , namely higher SG moments, g ( k ) gets smaller and the moments are more damped. Remark 3.12.
For multi-dimensional uncertainties ξ ∈ R N , we use the filter operator L ( u ( ξ )) = ( L ◦ · · · ◦ L N )( u ( ξ )) , where L n ( u ( ξ )) = ∂ ξ n (cid:0) (1 − ( ξ n ) (cid:1) ∂ ξ n u ( ξ ) . Finding a good choice of the filter strength λ , which sufficiently damps oscillations while preserving thesolution structure, is challenging since the optimal filter strength is problem dependent and a parameterstudy has to be performed.Within the derivation of the L filter, the use of the filter function has been motivated by an optimizationproblem of the form (3.14). Another approach is to replace the solution ansatz (3.13) by (cid:98) u K Ξ ( t, x, y, ξ ) := N Ξ (cid:88) l =1 K Ξ (cid:88) k =0 g (cid:18) kK Ξ (cid:19) λ ∆ t u k,l ( t, x, y ) φ k,l ( ξ ) χ l ( ξ ) , (3.16)Now, the function g : R + → R + plays the role of a filter function and its choice crucially affects the resultingapproximation quality. Let g fulfill the following properties:1. g (0) = 1,2. g ( (cid:96) ) (0) = 0 , for (cid:96) = 1 , . . . , p − g ( (cid:96) ) (1) = 0 , for (cid:96) = 0 , . . . , p − v ( s ) via v ( s ) ≈ (cid:98) v K Ξ ( s ) := K Ξ (cid:88) k =0 g (cid:18) kK Ξ (cid:19) a k φ k ( s ) , (3.17)where φ k ( s ) are basis functions and a k the corresponding coefficients of the expansion with k = 0 , . . . , K Ξ ,the filtered approximation (3.17) fulfills 9. If v has p continuous derivatives, then (cid:12)(cid:12) v ( s ) − (cid:98) v K Ξ ( s ) (cid:12)(cid:12) ≤ M · K Ξ1 / − p .
2. If v has a jump discontinuity at c ∗ ∈ [ − , (cid:12)(cid:12) v ( s ) − (cid:98) v K Ξ ( s ) (cid:12)(cid:12) ≤ M · d ( s ) − p K Ξ1 − p , where d ( s ) is the distance to c ∗ and M > p . The second inequality shows that even for discontinuous solutions, spectral accuracy is still possible if s is not too close to the shock position. One choice for the filter function in (3.16) is the exponential filterproposed in [21] which reads g (cid:18) kK Ξ (cid:19) = exp (cid:18) c (cid:18) kK Ξ (cid:19) α (cid:19) , (3.18)for k = 0 , . . . , K Ξ . Here, we use c = log( ε M ) where ε M is the machine accuracy. The parameter α is aninteger called the filter order. A more detailed discussion of the exponential filter can for example be foundin [8, Chapter 18.20]. We summarize the whole method for the filtered hyperbolicity-preserving stochastic Galerkin scheme in thefollowing algorithm, whereas we denote the solution operator as the right hand side of (3.9) by L ( n ) h ( u ( n ) k,i,j,l ) := u ( n ) k,i,j,l − ∆ t ∆ x (cid:90) Ξ l (cid:16) ˆ f ( u ( n ) i,j , u ( n ) i +1 ,j , ξ ) − ˆ f ( u ( n ) i − ,j , u ( n ) i,j , ξ ) (cid:17) φ k,l f Ξ l d ξ − ∆ t ∆ y (cid:90) Ξ l (cid:16) ˆ f ( u ( n ) i,j , u ( n ) i,j +1 , ξ ) − ˆ f ( u ( n ) i,j − , u ( n ) i,j , ξ ) (cid:17) φ k,l f Ξ l d ξ, (3.19)for k = 0 , . . . , K Ξ , ( i, j ) ∈ N X and l = 1 , . . . , N Ξ . Algorithm 1
Filtered Hyperbolicity-Preserving Stochastic Galerkin Scheme u (0) ← (cid:16) (cid:82) X i,j (cid:82) Ξ l u (0) φ k,l f Ξ l d ξ d( x, y ) (cid:17) ( i,j ) ∈ N X ,l =1: N Ξ ,k =0: K Ξ for n = 0 , . . . , N T do Choose the filter strength λ and filter order α (cid:98) u ( n ) ← (cid:16) g ( · ) u ( n ) k,i,j,l (cid:17) ( i,j ) ∈ N X ,l =1: N Ξ ,k =0: K Ξ (3.15) or (3.18) (cid:98) u ( n ) ← ΛΠ θK Ξ (cid:0)(cid:98) u ( n ) (cid:1) (3.10) (cid:98) u ( n +1) ← L ( n ) h ( (cid:98) u ( n ) ) (3.9) end for In our numerical examples we will apply the exponential filter (3.18).10 . Multi-element Intrusive Polynomial Moment Method
The multi-element approach from Section 3.1 is known to be applied to the stochastic Galerkin schemesince a local basis in the random variable yields improved solution approximations, even when using smallpolynomial degrees in each local element. Therefore, it is straightforward to adapt the strategy to the IPMmethod. In addition to that, such a local basis seems to be an ideal choice for IPM, since it decouples the IPMoptimization problems, enhancing the use of parallel implementations while heavily reducing computationalcosts.Given a strictly convex entropy h : R d → R of the uncertain system of hyperbolic conservation laws (2.1),the IPM system from [46] is given by ∂ t (cid:90) Ξ u φ ˜ k f Ξ d ξ + ∇ x · (cid:90) Ξ F (cid:18) ( ∇ u h ) − (cid:16) K Ξ (cid:88) k =0 λ k φ k (cid:17) , ξ (cid:19) φ ˜ k f Ξ d ξ = , (4.1)for ˜ k = 0 , . . . , K Ξ and where Λ := K Ξ (cid:88) k =0 λ k φ k is the entropic variable of the system, defining the solution u of (2.1) via the one-to-one map u = ( ∇ u h ) − ( Λ ).Performing the gPC expansion on the entropic variable guarantees hyperbolicity of the IPM system (4.1).For more information and the derivation of the system (4.1) we refer to [46]. The hyperbolicity of the IPMsystem is for example shown in [45, Chapter 4.2.3].We now follow the theory of [46] in order to obtain the IPM system (4.1) for the multi-element approachdescribed in Section 3.1. Thus, we divide the random domain Ξ into multi-elements Ξ l with l = 1 , . . . , , N Ξ and define a local basis φ k,l according to (3.2).We consider the ME system (3.6) and derive the corresponding IPM closure, i.e. we minimize the functional L ( u , λ ) = (cid:90) Ξ h ( u ) f Ξ d ξ + N Ξ (cid:88) l =1 K Ξ (cid:88) k =0 λ k,l (cid:16) u k,l − (cid:90) Ξ l u φ k,l f Ξ l d ξ (cid:17) , where we now have Lagrange multipliers λ ,l , . . . , λ K Ξ ,l ∈ R d for every multi-element Ξ l , l = 1 , . . . , N Ξ . By λ we denote these multipliers in vectorized form. In order to compute the exact minimizer of this functional,we calculate the Gˆateaux derivative with respect to u in an arbitrary direction v d u L ( u , λ )[ v ] := ddt L ( u + t · v , λ ) | t =0 . We obtain d u L ( u , λ )[ v ] = (cid:90) Ξ ∇ u h ( u ) · v f Ξ d ξ − N Ξ (cid:88) l =1 K Ξ (cid:88) k =0 (cid:90) Ξ l λ k,l φ k,l · v f Ξ l d ξ, which should be zero for any direction v , yielding ∇ u h ( u ) = N Ξ (cid:88) l =1 K Ξ (cid:88) k =0 λ k,l φ k,l . If we now define the entropic variable as Λ = N Ξ (cid:88) l =1 K Ξ (cid:88) k =0 λ k,l φ k,l ,
11e again have u = ( ∇ u h ) − ( Λ ). Plugging this ansatz into L and differentiating with respect to λ ˜ k, ˜ l for˜ k = 0 , . . . , K Ξ and ˜ l = 1 , . . . , N Ξ givesd λ ˜ k, ˜ l L ( u , λ ) = (cid:90) Ξ ˜ l ∇ u h ( u ( Λ )) ∇ Λ u ( Λ ) φ ˜ k, ˜ l f Ξ ˜ l d ξ + u ˜ k, ˜ l − (cid:90) Ξ ˜ l u ( Λ ) φ ˜ k, ˜ l f Ξ ˜ l d ξ − N Ξ (cid:88) l =1 K Ξ (cid:88) k =0 λ k,l (cid:90) Ξ ˜ l ∇ Λ u ( Λ ) φ k,l φ ˜ k, ˜ l f Ξ ˜ l d ξ = (cid:90) Ξ ˜ l Λ ∇ Λ u ( Λ ) φ ˜ k, ˜ l f Ξ ˜ l d ξ + u ˜ k, ˜ l − (cid:90) Ξ ˜ l u ( Λ ) φ ˜ k, ˜ l f Ξ ˜ l d ξ − (cid:90) Ξ ˜ l ∇ Λ u ( Λ ) Λ φ ˜ k, ˜ l f Ξ ˜ l d ξ = u ˜ k, ˜ l − (cid:90) Ξ ˜ l u ( Λ ) φ ˜ k, ˜ l f Ξ ˜ l d ξ ! = . (4.2)In order to minimize the Lagrangian we need to find the root of this equation. Note that since φ k,l hassupport Ξ l , the integral term becomes (cid:90) Ξ ˜ l u ( Λ ) φ ˜ k, ˜ l f Ξ ˜ l d ξ = (cid:90) Ξ ˜ l u (cid:32) N Ξ (cid:88) l =1 K Ξ (cid:88) k =0 λ k,l φ k,l (cid:33) φ ˜ k, ˜ l f Ξ ˜ l d ξ = (cid:90) Ξ ˜ l u (cid:32) K Ξ (cid:88) k =1 λ k, ˜ l φ k, ˜ l (cid:33) φ ˜ k, ˜ l f Ξ ˜ l d ξ. Hence, we can solve (4.2) through u ˜ k, ˜ l − (cid:90) Ξ ˜ l ( ∇ u h ) − (cid:32) K Ξ (cid:88) k =0 λ k, ˜ l φ k, ˜ l (cid:33) φ ˜ k, ˜ l f Ξ ˜ l d ξ = , (4.3)for ˜ k = 0 , . . . , K Ξ and ˜ l = 1 , . . . , N Ξ . In order to find the dual variables λ as the roots of (4.3) we employNewton’s method. Moreover, the ME-IPM system reads ∂ t (cid:90) Ξ l u φ ˜ k,l f Ξ l d ξ + ∇ x · (cid:90) Ξ l F (cid:0) Λ , ξ (cid:1) φ ˜ k,l f Ξ l d ξ = , (4.4)with ˜ k = 0 , . . . , K Ξ and l = 1 , . . . , N Ξ . Remark 4.1.
The IPM optimization problems decouple and one can solve the N Ξ problems per spatial cellindividually. Commonly, the multi-element approach allows using a smaller total degree of polynomials.Therefore, instead of solving one expensive optimization problem, we now solve N Ξ cheaper optimizationproblems in each cell, which can be distributed to different processors. We apply the numerical discretization from Section 3.2 yielding the following algorithm of the multi-elementintrusive polynomial moment method, whereas we now replace the solution operator (3.19) by L ( n ) h (cid:0) u ( n ) k,i,j,l , ¯ Λ ( n ) i,j,l (cid:1) := u ( n ) k,i,j,l − ∆ t ∆ x (cid:90) Ξ l (cid:16) ˆ f (cid:0) ( ∇ u h ) − ( ¯ Λ ( n ) i,j,l ) , ( ∇ u h ) − ( ¯ Λ ( n ) i +1 ,j,l ) , ξ (cid:1) − ˆ f (cid:0) ( ∇ u h ) − ( ¯ Λ ( n ) i − ,j,l ) , ( ∇ u h ) − ( ¯ Λ ( n ) i,j,l ) , ξ (cid:1)(cid:17) φ k,l f Ξ l d ξ − ∆ t ∆ y (cid:90) Ξ l (cid:16) ˆ f (cid:0) ( ∇ u h ) − ( ¯ Λ ( n ) i,j,l ) , ( ∇ u h ) − ( ¯ Λ ( n ) i,j +1 ,l ) , ξ (cid:1) − ˆ f (cid:0) ( ∇ u h ) − ( ¯ Λ ( n ) i,j − ,l ) , ( ∇ u h ) − ( ¯ Λ ( n ) i,j,l ) , ξ (cid:1)(cid:17) φ k,l f Ξ l d ξ, (4.5)12or k = 0 , . . . , K Ξ , ( i, j ) ∈ N X and l = 1 , . . . , N Ξ , defining one forward Euler time step of the FVM forthe ME-IPM system (4.4). Here, ¯ Λ ( n ) i,j,l describes the cell mean of the entropic variable in cell X i,j andmulti-element Ξ l at t n . Algorithm 2
Multi-Element Intrusive Polynomial Moment Method (ME-IPM) ¯ Λ (0) ← (cid:16) (cid:80) K Ξ k =0 (cid:82) X i,j (cid:82) Ξ l h (cid:48) (cid:0) u (0) (cid:1) φ k,l f Ξ l d ξ d( x, y ) φ k,l (cid:17) ( i,j ) ∈ N X ,l =1: N Ξ for n = 0 , . . . , N T do u ( n ) ← (cid:16) (cid:82) X i,j (cid:82) Ξ l ( ∇ u h ) − (cid:0) ¯ Λ ( n ) i,j,l (cid:1) φ k d ξ d( x, y ) (cid:17) ( i,j ) ∈ N X ,l =1: N Ξ ,k : K Ξ u ( n +1) ← L ( n ) h ( u ( n ) , ¯ Λ ( n ) ) (4.5) ¯ Λ ( n +1) ← arg min Λ L ( u ( n ) , ¯ Λ ( n ) ) (4.2) end forRemark 4.2. In addition to the described strategies to reduce computational costs, the steady-state test casein Section 5.2 makes use of two acceleration techniques taken from [31]. • Note that the minimization of the Lagrangian is usually performed using Newton iterations. Theiteration is stopped when the gradient of the Lagrangian is sufficiently close to zero. When solvingsteady-state problems, the FVM steps can be interpreted as a fixed point iteration, i.e. the time loopis iterated until a steady-state is reached. The idea of one-shot IPM (osIPM) is to perform only asingle Newton iteration to minimize the Lagrangian in every iteration of the time loop. In this case,we iterate both, the moments as well as the dual variables to their steady-state solution simultaneously.It can be shown that this iteration converges locally [31]. The idea follows the concept of the one-shotapproach used in PDE constraint optimization [18]. • We make use of adaptivity in stochastic space, which is one core advantage of intrusive methods[24, 31, 36, 54]. To ensure that a large number of time iterations is performed on a low refinementlevel (i.e. on a small number of moments), the maximal truncation order is only allowed to increaseonce the approximation is sufficiently close to the steady-state solution. Thereby, a large amountof time iterations use the lowest refinement level (i.e. a small number of moments), which can beunderstood as a preconditioning step. Keeping a low refinement level for a large amount of time stepsis called refinement retardation.
Remark 4.3.
A common issue of minimal entropy methods such as IPM is the loss of hyperbolicity. Thismeans that moment vectors generated by the numerical method can leave the hyperbolicity set, i.e. they donot belong to an admissible solution. The variable map ensures hyperbolicity as discussed in [26].
5. Numerical Results
In the following section we apply the oscillation mitigation intrusive UQ schemes from Section 3 and Section 4to the one- and two-dimensional Euler equations. In particular, we evaluate and compare the stochasticGalerkin and intrusive polynomial moment method with and without the multi-element ansatz in terms oftheir relative error and the computational costs.The code that we use to obtain the proceeding numerical results is openly available at [30] to allow repro-ducibility. 13 .1. Sod’s shock tube
To demonstrate the behavior of the different proposed methods, we first investigate the uncertain Sod’s shocktube test case as proposed in [46]. Here, the random system of conservation laws are the one-dimensionalEuler equations (3.12), where the heat capacity ratio γ is chosen to be 1 .
4. Initially, the gas is in a shockstate with random position x interface ( ξ ) = x + σξ , where ξ is uniformly distributed in the interval [ − , ρ IC = (cid:40) ρ L if x < x interface ( ξ ) ρ R else , ( ρv ) IC = 0 , ( ρe ) IC = (cid:40) ρ L e L if x < x interface ( ξ ) ρ R e R else . The IPM method chooses the entropy h ( ρ, ρv, ρe ) = − ρ ln (cid:18) ρ − γ (cid:18) ρe − ( ρv ) ρ (cid:19)(cid:19) . Further parameter values are X = [ x L , x R ] = [0 ,
1] range of spatial domain T = 0 .
14 end time x = 0 . , σ = 0 .
05 interface position parameters ρ L = 1 . , e L = 2 . , ρ R = 0 . , e R = 0 .
25 initial states
Table 1: Parameter values for Sod’s shock tube.
Moreover, we divide the spatial domain into N X = 2000 cells. In a first simulation, we compare the expectedgas density as well as its variance when employing stochastic Galerkin with and without the multi-elementansatz. To guarantee an admissible solution, both methods require the hyperbolicity-preserving limiterdiscussed in Section 3.3. The classical hSG scheme uses polynomials up to a degree of 14, which we denoteby hSG . For multi-element hSG we decompose the stochastic domain into three multi-elements, each ofthem with polynomials up to degree 4. We denote this method by ME-hSG , . Furthermore, a Gauss-Legendre quadrature is chosen. Here, the chosen number of points is 30 for the classical hSG method and10 points per multi-element when using ME-hSG.The results are shown in Figure 1. From left to right, we can see the expected value and variance of the threemain characteristics of this test case, namely a rarefaction wave, a contact discontinuity and a shock. Sincethe number of unknowns per spatial cell is identical for both hSG and ME-hSG, the computational costs arethe same leading to similar runtimes (24.675 seconds for hSG and 23.761 seconds for ME-hSG). Expectedvalue and variance of both methods look similar. In particular, the solution emits a step-like profile of theexpected value and an oscillatory variance at the shock. In this case, the hSG and ME-hSG results almostcoincide.When applying the IPM method with and without multi-elements we observe a clear improvement of theoverall runtime. We denote the IPM method using polynomials up to degree 14 by IPM and the multi-element IPM method with three multi-elements and gPC degree four by ME-IPM , , where we use the samequadrature as for hSG and ME-hSG. The numerical results are depicted in Figure 2.First, we notice that the multi-element ansatz heavily reduces the computational costs, since every cellrequires solving three decoupled optimization problems with five unknowns respectively. In contrast, theclassical IPM method with the same number of unknowns requires solving one optimization problem with 1514 .3 0.4 0.5 0.6 0.7 0.8x0.20.40.60.81.0 E[ ] reference hSG ME hSG reference hSG ME hSG
Figure 1: Comparison of hSG using 15 moments and ME-hSG for 3 multi-elements with 5 moments. The number of unknownsper spacial cell is identical, which yields the same runtime and similar results. referenceIPM ME IPM referenceIPM ME IPM
Figure 2: Comparison of IPM using 15 moments and ME-IPM for 3 multi-elements with 5 moments. Even though the numberof unknowns per spacial cell is identical, the runtime of the multi-element IPM method is reduced by a factor of 7 .
72 whileyielding a similar result. unknowns. This reduction in numerical costs is reflected by the runtime which decreases by a factor of 7 . λ = 2 andthe order α = 10. We observe that the filter significantly improves the expected value and the varianceapproximation at the shock. Unfortunately, the filter slightly dampens the variance at the rarefaction wave.Overall, the filter yields a satisfactory solution approximation without increasing the computational time.15 .3 0.4 0.5 0.6 0.7 0.8x0.20.40.60.81.0 E[ ] referenceME hSG ME fhSG referenceME hSG
ME fhSG
Figure 3: Comparison of ME-hSG and ME-fhSG using 3 multi-elements with 5 moments. While the filter slightly dampensthe variance at the contact discontinuity and the rarefaction wave, the shock approximation shows good agreement with thereference solution. referenceME hSG
ME fhSG referenceME hSG
ME fhSG
Figure 4: Results for ME-hSG and ME-fhSG using 3 multi-elements with 5 moments. Zoom on shock.
In the following we consider the uncertain NACA test case from [31]. The test case investigates the effectsof an uncertain angle of attack φ ∼ U (0 . , .
75) for a NACA0012 airfoil with a length of one meter. Thegoverning equations are the stochastic Euler equations in two dimensions, which read ∂ t ρρv ρv ρe + ∂ x ρv ρv + pρv v v ( ρe + p ) + ∂ y ρv ρv v ρv + pv ( ρe + p ) = . (5.1)16 closure for the pressure p in two dimensions is given by p = ( γ − ρ (cid:18) e −
12 ( v + v ) (cid:19) . The hyperbolicity set from Definition 3.3 is therefore given by R = u = ρρv ρv ρe (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ρ > , p = ( γ − ρ (cid:18) e −
12 ( v + v ) (cid:19) > . We again choose the heat capacity ratio γ as 1 .
4. At the airfoil surface, the Euler slip condition v T n = 0 isused, where n denotes the surface normal. At a distance of 40 meters from the airfoil, one assumes Dirichletboundary conditions, i.e. the airfoil does not influence the flow at this boundary. Here, the constant flowfield is given by Mach number M a = 0 .
8, pressure p = 101 325 Pa and a temperature of 273 .
15 K. With thespecific gas constant R = 287 . Jkg · K and a given angle of attack φ the conserved variables at the far fieldboundary can be determined uniquely, since ρ = pR · T ,ρv = ρ · Ma · (cid:112) γRT cos φ,ρv = ρ · Ma · (cid:112) γRT sin φ,ρe = pγ − ρ ( v + v ) . The initial condition is chosen to be equal to the far field boundary values. In this case, the slip conditionat the airfoil will correct the flow solution until a steady-state flow is reached after a certain time. Theiteration in time is performed until the time residual fulfills the stopping criterion (cid:88) ( i,j ) ∈ N X N Ξ (cid:88) l =1 ∆ x ∆ y (cid:107) u ( n ) i,j,l − u ( n − i,j,l (cid:107) ≤ ε. (5.2)In all steady-state calculations, we choose ε = 5 · − .Unlike the Shock tube, this test case does not have an analytical solution. Therefore, we compute a referencesolution by stochastic collocation with 100 Gauss-Legendre quadrature points. The resulting expected valueand variance of the density as well as the implemented mesh are depicted in Figure 5. The mesh is composedof 22361 triangular elements and is used for all computations of the NACA test case. Remark 5.1.
Within the numerical implementation of this test case, we divide the two-dimensional spatialdomain into a triangular mesh. The theory of Section 3.2 for uniform grids can easily be applied to such adiscretization.
In the following, we aim at approximating this solution with our oscillation mitigation intrusive UQ methods.To quantify the efficiency of the different approaches, we investigate the relative error of the density’sexpected value and variance over time. We define the L -error of a discrete quantity e ∆ ∈ R N x × N y (cid:107) e ∆ (cid:107) ∆ := (cid:115) (cid:88) ( i,j ) ∈ N X ∆ x ∆ y e i,j . Then, if we denote the solution obtained with the numerical method by u ∆ and the reference solution by (cid:98) u ∆ , the relative error of the expected value and variance is (cid:107) E[ (cid:98) u ∆ ] − E[ u ∆ ] (cid:107) ∆ (cid:107) E[ (cid:98) u ∆ ] (cid:107) ∆ and (cid:107) Var[ (cid:98) u ∆ ] − Var[ u ∆ ] (cid:107) ∆ (cid:107) Var[ (cid:98) u ∆ ] (cid:107) ∆ . (5.3)17 igure 5: Left: Reference solution E[ ρ ]. Center: Var[ ρ ]. Right: Mesh around the airfoil. The mesh is a circle around the airfoil,consisting of 22361 triangular cells. Following [31], we record this error inside a box of one meter height and 1.1 meters length around the airfoil.This strategy allows us to exclude fluctuations in the far field solution. The resulting efficiency curves aredepicted in Figure 6. Here, we plot the relative errors according to (5.3) over the computational time.All methods are converged to a time residual of ε = 5 · − , cf. (5.2). For most methods, a satisfactorysolution is however reached at an earlier stage when the efficiency curve in Figure 6 is saturated. In thiscase, the error of the steady-state approximation reaches the same error level as the approximation error ofthe uncertainty. Error E[ ] hSG osIPM ME-hSG osME-IPM Error Var[ ] hSG osIPM ME-hSG osME-IPM
Figure 6: Relative errors (5.3) computed over runtime. Computations are carried out using 5 MPI threads. As for theone-dimensional shock tube, the multi-element ansatz significantly reduces the overall runtime when using the IPM method.Multi-element IPM beats multi-element hSG, however the fastest method yielding the smallest error is hSG.
All simulations have been computed in parallel on 5 processors. For details on the parallelization of the code,see [31]. The stochastic Galerkin calculations again require the hyperbolicity limiter presented Section 3.3.First, we compare the results obtained with IPM. Here, we make use of one-shot IPM (osIPM) to reduce18umerical costs. More information on this method can be found in [31]. The number of unknowns perspatial cell for multi-element IPM and IPM is chosen to be 9: While ME-IPM employs 3 multi-elementswith gPC degree two, the classical IPM method has gPC degree 8. We use a Clenshaw-Curtis quadratureof level 4 (17 points) when not using the multi-element ansatz. The multi-element methods use a level 2 (5points) Clenshaw-Curtis rule per element. As expected, the multi-element ansatz again heavily decreasesthe numerical costs, while yielding a similar error level. Here, the error obtained with the standard IPMmethod is slightly better compared to the multi-element error. Note that ME-IPM beats ME-hSG in termsof efficiency. In this test case the most efficient method is the hyperbolicity-preserving stochastic Galerkinapproach.In order to further decrease the error obtained with the multi-element methods, we make use of adaptivityand refinement retardation, cf. Remark 4.2. Here, we fix the choice of 3 elements, however let the numberof polynomials in every element depend on the solution’s smoothness. Furthermore, all quadratures areapproximated by a Clenshaw-Curtis quadrature rule, which adapts to the moment order. The differentmoment orders and the corresponding quadrature points for each refinement level are given by table Table 2.We choose the refinement barriers as δ − = 2 · − and δ + = 2 · − . For more information on Table 2 andrefinement barriers see [31]. refinement level 0 1 2 3 4moment order 2 3 4 5 6quadrature points 5 9 9 9 17 Table 2: Refinement levels with their moment orders and quadrature points for adaptive refinements in ME-hSG and ME-IPM.
We converge the solution on refinement level zero to a residual of ε = 4 · − and then let the refinementlevel allow to increase up to level 4 until ε = 5 · − . The resulting efficiency is depicted in Figure 7. In thiscase, the ME stochastic Galerkin method is more efficient than ME-IPM. This is most likely caused by theincreased number of moments per multi-element, which increases the costs of solving the IPM optimizationproblem. Thus, the ME-IPM method is most efficient on low polynomial degrees. A refinement strategy,which overcomes this issue would be to increase or decrease the number of multi-elements instead of thenumber of polynomials. We however leave this idea to future work. Error E[ ] readosME-IPM readME-hSG Error Var[ ] readosME-IPM readME-hSG
Figure 7: Relative errors (5.3) computed over runtime. Computations are carried out using 5 MPI threads. We use refinementretardation (re) and adaptivity (ad) to accelerate the computation while achieving an improved error level.
Figure 8 shows the approximated expected value and variance by refinement retardation adaptive ME-hSG(readME-hSG) as well as the corresponding refinement level at the steady-state solution. The results ofreadosME-IPM look similar and are therefore left out at this point. It can be seen that the approximation19oes not show oscillatory behaviour, which stems from the sufficiently high number of unknowns chosenby the adaptive scheme on top of the airfoil. The approximation agrees well with the reference solution inFigure 5.
Figure 8: Expected density (left) and corresponding variance (right) computed with ME-hSG using refinement retardation andadaptivity. Right: Refinement level at final solution.
We leave the use of filters for steady-state problems to future work. Here, the main issue is that the filteris applied in every pseudo-time step, which commonly leads to a heavy modification of the solution.
Lastly, we investigate a nozzle with an uncertain shock. Again, the test case relies on the two-dimensionalEuler equations (5.1) inside a nozzle geometry (see Figure 9). The nozzle is composed of a chamber onthe left, a throat in the center and the main nozzle part on the right. The shock position is uniformlydistributed within the throat. The shock states are identical to the one-dimensional shock tube experimentfrom Section 5.1, namely ρ IC = (cid:40) ρ L if x < x interface ( ξ ) ρ R else , ( ρv ) IC = 0 , ( ρv ) IC = 0 , ( ρe ) IC = (cid:40) ρ L e L if x < x interface ( ξ ) ρ R e R else , with x interface ( ξ ) = x + σξ , ξ ∼ U ( − , x to be the x-coordinate of the center of the throat.The throat ranges from x − . x + 0 .
5, i.e. we set σ = 0 . t = 4. A reference solution,which has been computed by stochastic collocation with a 100 point Gauss-Legendre quadrature, is shownin Figure 9. All computations are carried out on a mesh with N X = 225758 triangular cells. Note thatsimilar to Sod’s shock tube we observe structures stemming from the rarefaction wave (left), the contactdiscontinuity (center) and the shock (right). 20 igure 9: Reference solution for density. Top: Expected value. Bottom: Variance. We first compare the different methods without multi-elements. Here, we use gPC polynomials up to degree5 as well as a 20 point Gauss-Legendre quadrature rule. We run simulations with hSG, fhSG and IPM. Thefiltered hSG method is implemented for the exponential filter with order α = 7 and strength λ = 0 . igure 10: Expected density with different methods. From top to bottom: hSG, fhSG, IPM. As before, we apply the exponential filter with order α = 7 and strength λ = 0 .
1. When using multi-elements, the filter is able to smear out oscillations, leading to an improved approximation of the expectedvalue and the variance. However, the effectiveness of the filter appears to be increased when applying it tothe standard hSG method. Furthermore, note that all multi-element methods yield a better resolution ofthe variance inside the throat region. 22 igure 11: Variance of the density with different methods. From top to bottom: hSG, fhSG, IPM. igure 12: Expected density with different methods. From top to bottom: ME-hSG, ME-fhSG, ME-IPM. igure 13: Variance of the density with different methods. From top to bottom: ME-hSG, ME-fhSG, ME-IPM. . Conclusions and Outlook In this work, we combined different methods to mitigate oscillations and spurious artifacts that arise fromintrusive UQ methods while maintaining hyperbolicity of the moment system. To mitigate oscillations, wemade use of filtering as well as a multi-element ansatz. Furthermore, hyperbolicity is ensured by usinga hyperbolicity limiter as well as the IPM method. First, we combined filtering with the hyperbolicity-preserving limiter and then extended the resulting method with the multi-element ansatz. Second, themulti-element ansatz has been applied to the IPM method. Both strategies dampen oscillations which leadsto an improved approximation of expected values and variances. In addition to that, the multi-elementapproach in combination with IPM allows heavily reducing the required runtime to achieve a certain errorlevel. When comparing the different methods, we observe that the optimal method choice is problemdependent. In the case of Sod’s shock tube, filtering in combination with the multi-element ansatz in hSGyields the best approximation. For the nozzle test case, applying the filter to the original hSG methodappears to be the method of choice as it closely captures the reference solution. A solely application ofthe (necessary) hyperbolicity limiter or the IPM closure to this test case resulted in oscillatory solutionapproximations. However, either the use of the multi-element ansatz or the filter significantly improves theresolution quality. The presented NACA test case underlines the efficiency of the multi-element approachin combination with IPM.The different methods are so far applied for first-order numerical schemes and it would be interesting to becombined with WENO reconstructions in order to obtain high-order approximations and to further reducethe impact of Gibbs phenomenon. In addition to that, the filters that we presented within this article can beused together with IPM in order to filter the coefficients of the entropic variable. This might yield a furtherreduction of oscillations within IPM. Since it is not clear how to achieve hyperbolicity of filtered moments,we leave this task to future work. Furthermore, the applicability of filters for steady problems should beinvestigated to allow the use of filters in steady state applications.
Acknowledgments
Funding by the Deutsche Forschungsgemeinschaft (DFG) within the RTG GrK 1932 “Stochastic Models forInnovations in the Engineering Science” is gratefully acknowledged. Jonas Kusch has been supported by theDFG under grant FR 2841/6-1.
References [1]
R. Abgrall , A simple, flexible and generic deterministic approach to uncertainty quantifications in non-linear problems:application to fluid flow problems , tech. rep., 2008.[2]
R. Abgrall and P. M. Congedo , A semi-intrusive deterministic approach to uncertainty quantification in non-linearfluid flow problems , Journal of Computational Physics, 235 (2013), pp. 828–845.[3]
R. Abgrall and S. Mishra , Uncertainty quantification for hyperbolic systems of conservation laws , in Handbook ofNumerical Analysis, vol. 18, Elsevier, 2017, pp. 507–544.[4]
B. K. Alpert , A Class of Bases in L2 for the Sparse Representation of Integral Operators , SIAM Journal on MathematicalAnalysis, 24 (1993), pp. 246–262.[5]
T. Barth , On the propagation of statistical model parameter uncertainty in CFD calculations , Theoretical and Compu-tational Fluid Dynamics, 26 (2012), pp. 435–457.[6] ,
Non-intrusive uncertainty propagation with error bounds for conservation laws containing discontinuities , inUncertainty quantification in computational fluid dynamics by H. Bijl, D. Lucor, S. Mishra, C. Schwab, Springer, 2013,pp. 1–57.[7]
J. P. Boyd , The erfc-log filter and the asymptotics of the Euler and Vandeven sequence accelerations , in Proceedings ofthe Third International Conference on Spectral and High Order Methods, Houston Math. J, 1996, pp. 267–276.[8]
J. P. Boyd , Chebyshev and Fourier spectral methods , Courier Corporation, 2001.[9]
Q. Y. Chen, D. Gottlieb, and J. S. Hesthaven , Uncertainty analysis for the steady-state flows in a dual throat nozzle ,Journal of Computational Physics, 204 (2005), pp. 378–398. B. J. Debusschere, H. N. Najm, P. P. P´ebay, O. M. Knio, R. G. Ghanem, and O. P. Le Maıtre , Numerical challengesin the use of polynomial chaos representations for stochastic processes , SIAM journal on scientific computing, 26 (2004),pp. 698–719.[11]
S. Deshpande , Kinetic theory based new upwind methods for inviscid compressible flows , in 24th Aerospace SciencesMeeting, 1986, p. 275.[12]
B. Despr´es, G. Po¨ette, and D. Lucor , Robust Uncertainty Propagation in Systems of Conservation Laws with theEntropy Closure Method , Springer International Publishing, 2013, pp. 105–149.[13]
J. D¨urrw¨achter, T. Kuhn, F. Meyer, L. Schlachter, and F. Schneider , A hyper-bolicity-preserving discontinu-ous stochastic Galerkin scheme for uncertain hyperbolic systems of equations , Journal of Computational and AppliedMathematics, 370 (2020), p. 112602.[14]
M. B. Giles , Multilevel Monte Carlo path simulation , Operations Research, 56 (2008), pp. 607–617.[15]
D. Gottlieb and D. Xiu , Galerkin method for wave equations with uncertain coefficients , Communications in Compu-tational Physics, 3 (2008), pp. 505–518.[16]
A. Harten, P. D. Lax, and B. v. Leer , On upstream differencing and godunov-type schemes for hyperbolic conservationlaws , SIAM review, 25 (1983), pp. 35–61.[17]
C. D. Hauck , High-order entropy-based closures for linear transport in slab geometry , Communications in MathematicalSciences, 9 (2011), pp. 187–205.[18]
S. B. Hazra, V. Schulz, J. Brezillon, and N. R. Gauger , Aerodynamic shape optimization using simultaneous pseudo-timestepping , Journal of Computational Physics, 204 (2005), pp. 46–64.[19]
S. Heinrich , Multilevel Monte Carlo methods , Large-scale scientific computing, Third international conference LSSC 2001,Sozopol, Bulgaria, 2170 (2001), pp. 58–67.[20]
J. S. Hesthaven, S. Gottlieb, and D. Gottlieb , Spectral methods for time-dependent problems , vol. 21, CambridgeUniversity Press, 2007.[21]
B. J. Hoskins , Representation of the earth topography using spherical harmonies , Monthly Weather Review, 108 (1980),pp. 111–115.[22]
S. Jin and L. Pareschi , eds.,
Uncertainty Quantification for Hyperbolic and Kinetic Equations , vol. 14 of SEMA SIMAISpringer Series, Springer International Publishing, Cham, 2017.[23]
O. Kolb , A Third Order Hierarchical Basis WENO Interpolation for Sparse Grids with Application to ConservationLaws with Uncertain Data , Journal of Scientific Computing, 74 (2018), pp. 1480–1503.[24]
I. Kr¨oker and C. Rohde , Finite volume schemes for hyperbolic balance laws with multiplicative noise , Applied NumericalMathematics, 62 (2012), pp. 441–456.[25]
J. Kusch , Uncertainty Quantification for Hyperbolic Equations , RWTH Aachen University, (2015), pp. 1–23.[26]
J. Kusch, G. W. Alldredge, and M. Frank , Maximum-principle-satisfying second-order Intrusive Polynomial Momentscheme , The SMAI journal of computational mathematics, 5 (2019), pp. 23–51.[27]
J. Kusch and M. Frank , Intrusive methods in uncertainty quantification and their connection to kinetic theory , Inter-national Journal of Advances in Engineering Sciences and Applied Mathematics, 10 (2018), pp. 54–69.[28] ,
An adaptive quadrature-based moment closure , International Journal of Advances in Engineering Sciences andApplied Mathematics, 11 (2019), pp. 174–186.[29]
J. Kusch, R. G. McClarren, and M. Frank , Filtered Stochastic Galerkin Methods For Hyperbolic Equations , Journalof Computational Physics, 403 (2020).[30]
J. Kusch, L. Schlachter, and J. Wolters , UQCreator testcases for ”Oscilla-tion Mitigation of Hyperbolicity-Preserving Intrusive Uncertainty Quantification Meth-ods for Systems of Conservation Laws” , 2020. https://git.scc.kit.edu/uqcreator/publication-oscillation-mitigation-of-intrusive-uncertainty-quantification-methods-for-hyperbolic-systems .[31]
J. Kusch, J. Wolters, and M. Frank , Intrusive acceleration strategies for uncertainty quantification for hyperbolicsystems of conservation laws , Journal of Computational Physics, (2020), p. 109698.[32]
C. D. Levermore , Moment closure hierarchies for kinetic theories , Journal of statistical Physics, 83 (1996), pp. 1021–1065.[33]
K. O. Lye , Multilevel Monte-Carlo for measure valued solutions , (2016).[34]
X. Ma and N. Zabaras , An adaptive hierarchical sparse grid collocation algorithm for the solution of stochastic differentialequations , Journal of Computational Physics, 228 (2009), pp. 3084–3113.[35]
R. G. McClarren and C. D. Hauck , Robust and accurate filtered spherical harmonics expansions for radiative transfer ,Journal of Computational Physics, 229 (2010), pp. 5597–5614.[36]
F. Meyer, C. Rohde, and J. Giesselmann , A posteriori error analysis for random scalar conservation laws using thestochastic galerkin method , IMA Journal of Numerical Analysis, 40 (2020), pp. 1094–1121.[37]
S. Mishra, N. H. Risebro, C. Schwab, and S. Tokareva , Numerical solution of scalar conservation laws with randomflux functions , SIAM-ASA Journal on Uncertainty Quantification, 4 (2016), pp. 552–591.[38]
S. Mishra and C. Schwab , Sparse tensor multi-level Monte Carlo finite volume methods for hyperbolic conservation lawswith random initial data , Mathematics of Computation, 81 (2012), pp. 1979–2018.[39]
S. Mishra, C. Schwab, and J. ˇSukys , Multi-level Monte Carlo finite volume methods for nonlinear systems of conser-vation laws in multi-dimensions , Journal of Computational Physics, 231 (2012), pp. 3365–3388.[40] ,
Multi-level Monte Carlo Finite Volume Methods for Uncertainty Quantification in Nonlinear Systems of BalanceLaws , Lecture Notes in Computational Science and Engineering, 92 (2013).[41]
F. Nobile, R. Tempone, and C. G. Webster , A Sparse Grid Stochastic Collocation Method for Partial DifferentialEquations with Random Input Data , SIAM Journal on Numerical Analysis, 46 (2008), pp. 2411–2442. B. Perthame , Boltzmann type schemes for gas dynamics and the entropy property , SIAM Journal on Numerical Analysis,27 (1990), pp. 1405–1421.[43]
B. Perthame , Second-order boltzmann schemes for compressible euler equations in one and two space dimensions , SIAMJournal on Numerical Analysis, 29 (1992), pp. 1–19.[44]
B. Perthame and C.-W. Shu , On positivity preserving finite volume schemes for Euler equations , Numerische Mathe-matik, 73 (1996), pp. 119–130.[45]
G. Po¨ette , Contribution to the mathematical and numerical analysis of uncertain systems of conservation laws and ofthe linear and nonlinear Boltzmann equation , PhD thesis, 2019.[46]
G. Po¨ette, B. Despr´es, and D. Lucor , Uncertainty quantification for systems of conservation laws , Journal of Com-putational Physics, 228 (2009), pp. 2443–2467.[47]
G. Po¨ette, B. Despr´es, and D. Lucor , Treatment of uncertain material interfaces in compressible flows , ComputerMethods in Applied Mechanics and Engineering, 200 (2011), pp. 284–308.[48]
D. Radice, E. Abdikamalov, L. Rezzolla, and C. D. Ott , A new spherical harmonics scheme for multi-dimensionalradiation transport I. Static matter configurations , Journal of Computational Physics, 242 (2013), pp. 648–669.[49]
C. Schillings and A. M. Stuart , Analysis of the ensemble Kalman filter for inverse problems , SIAM Journal onNumerical Analysis, 55 (2017), pp. 1264–1290.[50]
L. Schlachter and F. Schneider , A hyperbolicity-preserving stochastic Galerkin approximation for uncertain hyperbolicsystems of equations , Journal of Computational Physics, 375 (2018), pp. 80–98.[51]
L. Schlachter, F. Schneider, and O. Kolb , Weighted Essentially Non-Oscillatory stochastic Galerkin approximationfor hyperbolic conservation laws , Journal of Computational Physics, 419 (2020), p. 109663.[52]
F. Schneider, J. Kall, and G. W. Alldredge , A realizability-preserving high-order kinetic scheme using WENOreconstruction for entropy-based moment closures of linear kinetic equations in slab geometry , Kinetic and Related Models,9 (2015), pp. 193–215.[53]
J. Tryoen and A. Ern , Adaptive Anisotropic Spectral Stochastic Methods , (2010).[54]
J. Tryoen, O. L. Maıtre, and A. Ern , Adaptive anisotropic spectral stochastic methods for uncertain scalar conservationlaws , SIAM Journal on Scientific Computing, 34 (2012), pp. A2459–A2481.[55]
L. M. M. van den Bos, B. Sanderse, W. A. A. M. Bierbooms, and G. J. W. van Bussel , Bayesian model calibrationwith interpolating polynomials based on adaptively weighted Leja nodes , Communications in Computational Physics, 27(2018), pp. 33–69.[56]
X. Wan and G. E. Karniadakis , An adaptive multi-element generalized polynomial chaos method for stochastic differ-ential equations , Journal of Computational Physics, 209 (2005), pp. 617–642.[57] ,
Long-term behavior of polynomial chaos in stochastic flow simulations , Computer Methods in Applied Mechanicsand Engineering, 195 (2006), pp. 5582–5596.[58] ,
Multi-Element Generalized Polynomial Chaos for Arbitrary Probability Measures , SIAM J. Sci. Comput., 28(2006), pp. 901–928.[59] ,
Error control in multi-element generalized polynomial chaos method for elliptic problems with random coefficients ,Communications in Computational Physics, 5 (2009), pp. 793–820.[60]
N. Wiener , The homogeneous chaos. , Amer. J. Math, 60 (1938), pp. 897–936.[61]
J. A. Witteveen, A. Loeven, and H. Bijl , An adaptive Stochastic Finite Elements approach based on Newton-Cotesquadrature in simplex elements , Computers and Fluids, 38 (2009), pp. 1270–1288.[62]
D. Xiu and J. S. Hesthaven , High-order collocation methods for differential equations with random inputs , 27 (2005),pp. 1118–1139.[63]
D. Xiu and G. E. Karniadakis , The Wiener-Askey polynomial chaos for stochastic differential equations , SIAM Journalon Scientific Computing, 24 (2003), pp. 619–644.[64]
L. Yan and T. Zhou , Adaptive multi-fidelity polynomial chaos approach to Bayesian inference in inverse problems ,Journal of Computational Physics, 381 (2019), pp. 110–128.[65]
X. Zhang and C.-W. Shu , On positivity preserving high order discontinuous Galerkin schemes for compressible Eulerequations on rectangular meshes , Journal of Computational Physics, 229 (2010), pp. 8918–8934., Journal of Computational Physics, 229 (2010), pp. 8918–8934.