The use of spectral projection-based methods for simulation of a stochastic system with discontinuous solution exhibits the Gibbs phenomenon, which is characterized by oscillations near discontinuities. This paper investigates a dynamic bi-orthogonality-based approach with appropriate postprocessing for mitigating the effects of the Gibbs phenomenon. The proposed approach uses spectral decomposition of the spatial and stochastic fields in appropriate orthogonal bases, while the dynamic orthogonality (DO) condition is used to derive the resultant closed-form evolution equations. The orthogonal decomposition of the spatial field is exploited to propose a Gegenbauer reprojection-based postprocessing approach, where the orthogonal bases in spatial dimension are reprojected on the Gegenbauer polynomials in the domain of analyticity. The resultant spectral expansion in Gegenbauer series is shown to mitigate the Gibbs phenomenon. Efficacy of the proposed method is demonstrated for simulation of a one-dimensional stochastic Burgers equation and stochastic quasi-one-dimensional flow through a convergent-divergent nozzle.
Continuous advancements in digital technologies have resulted in ubiquitous use of computer simulators for investigation of many large-scale systems. The simulator predictions are often uncertain due to unknown/poorly known physics, model parameters, initial and boundary conditions. Need and importance of uncertainty quantification in simulation predictions has already been emphasized by researchers in varied fields [1–4]. Subsequently, the research community has invested significant efforts to develop various uncertainty quantification methodologies (see Ref.  and references therein.). In essence, uncertainty quantification is an estimate of joint probability distribution of system responses conditional on probabilistic specification of uncertainties in the simulator model, parameters, initial and boundary conditions, etc.
where is a differential operator.
Owing to the generality and ease of implementation, Monte Carlo methods are widely used for solution of Eq. (1). The method involves simulation of the system at large number of conditions sampled from input probability distribution. The efficiency of the Monte Carlo method can be improved by using variance reduction techniques like importance sampling  or sequential Monte Carlo methods , among others. Note that the Monte Carlo methods exhibit the accuracy proportional to , where Ns is the number of samples collected. Thus, the Monte Carlo method typically requires a large number of samples for acceptable accuracy . For the simulation of large-scale systems, the computational cost may render application of Monte Carlo methods for uncertainty quantification intractable.
Note that choice of the orthogonal basis is made based on the Hilbert spaces and . A Karhunen–Loeve expansion  uses the tensor product (5), where is spectrally represented in terms of orthogonal basis in with temporally varying expansion coefficients in . A generalized polynomial chaos (gPC) expansion  is an example of tensor product (7) with orthogonal basis in and the associated expansion coefficients in . This paper investigates a bi-orthogonal method, which uses the gPC basis in , while the dynamically orthogonal eigenfunction basis is used in .
The present work is motivated by the need to quantify uncertainty in simulation of systems with discontinuous solutions. An example of such a discontinuity is shocks that evolve during transonic simulations. In particular, the paper deals with stochastic systems defined using the hyperbolic conservation laws that invariably result in a discontinuous solution irrespective of the initial conditions. The simulation of such system has attracted significant interest by the research community and the literature is rich with methods for resolution of discontinuities in deterministic simulations (see Ref.  and references therein for a review of the methods). However, the problem of the development of numerical methods for solution of stochastic hyperbolic partial differential equations is not yet resolved and is of current interest to the research community .
Considering the cost of the Monte Carlo methods, the present paper focuses on SSP methods for the solution of SPDE in Eq. (1). The use of spectral methods for hyperbolic SPDEs is known to pose the following significant numerical challenges 
Numerical solutions using finite difference schemes lead to spurious oscillations in the expansion coefficients and/or orthogonal basis.
The spectral expansion of the solution of hyperbolic SPDEs lead to Gibbs phenomenon as discontinuities develop in solution . The Gibbs phenomenon for application of SSP methods to hyperbolic SPDEs with discontinuous solution is characterized by
Slow convergence of the spectral approximation of the solution at points away from the discontinuity. The convergence is , where N is the number of expansion terms retained in the spectral approximation.
oscillations are observed in the solution near discontinuities that do not decrease with increasing N.
The existence and resolution of Gibbs phenomenon in spectral expansion of a deterministic function using a finite number of orthogonal bases is extensively investigated in the literature [19–21]. The spectral expansion of a function using global orthogonal basis is contaminated by the presence of discontinuities. However, the expansion coefficients contain enough information to recover the function with high accuracy using appropriate postprocessing. Gibbs phenomenon can be completely resolved by reprojecting the partial sum on Gibbs complementary basis . Gottlieb et al. [20–23] have shown that exponentially convergent approximation can be obtained at point values in a subinterval of analyticity by reprojecting the partial sum on orthogonal Gegenbauer polynomial basis. The postprocessing method can be extended to stochastic functions, provided, the spectral expansion is available in terms of orthogonal basis in spatial dimensions.
where is the mean, are eigenfunctions which form a complete orthogonal basis in at a given time step t, while are independent zero-mean random variables. Sapsis and Lermusiaux  have proposed the dynamic orthogonality (DO) condition, which constrains to be orthogonal to , to ensure orthogonality of eigenfunctions at all time steps.
where are orthogonal polynomials. The DO condition is used to derive closed-form solution of evolution equations for , and . However, the resultant approximation of in Eq. (10) using the solution of and exhibits the Gibbs phenomenon. A Gegenbauer reprojection-based postprocessing method is proposed to mitigate the Gibbs phenomenon. During the postprocessing step, eigenfunctions are reprojected on the Gegenbauer basis. Note that the Bi-orthogonal expansion is explored by Cheng et al. , where the independence of expansion coefficients is used to obtain Bi-orthogonality. Choi et al.  combine the gPC expansion with the DO method, where the gPC method is used during the initial time evolution, while the DO expansion is used subsequently when the uncertainty in the solution is high enough. Key contribution of the present paper is a Gegenbauer reprojection method for mitigation of the Gibbs Phenomenon.
There have been recent research efforts that focused on addressing the issue of Gibbs phenomenon in gPC settings. Wan and Karniadakis  have proposed a multi-element gPC (ME-gPC) method to mitigate the effect of Gibbs phenomenon (also see Lin et al. ). The ME-gPC method uses domain decomposition of the stochastic space with locally defined gPC basis, while the forward problem is solved for each subdomain. The ME-gPC method uses localized expansion in the subdomains to mitigate the Gibbs phenomenon. However, the computational cost of the ME-gPC method rises quickly as the number of subdomains increases. Poette et al.  have proposed an entropy-based intrusive polynomial moment method for uncertainty propagation in nonlinear systems with discontinuous solutions. This method exhibited better computationally efficient than the ME-gPC, but the requirement for minimization of the entropy at each time step results in the computational cost higher than the gPC method. The method reformulates the SPDE in terms of appropriately selected entropic variables that are bijection of the uncertain parameters, while the gPC expansion of the entropic variables is defined through the Galerkin projection of the bijection. The method is implemented in two steps: in the first step, the resultant stochastic Galerkin system is discretized using a finite volume approach, and upwinded row scheme is used for numerical solution. In the second step, gPC expansion of the entropic variables is calculated so that the entropy of the system is minimized. The method is demonstrated for solution of the stochastic Burgers equation. The method is computationally more efficient than the ME-gPC; however, the requirement for minimization of the entropy at each time step results in a computational cost higher than the gPC method. Tryoen et al.  use a Roe-type solver for intrusive Galerkin projection of uncertain hyperbolic systems. The solver uses the Dubois and Mehlman entropy corrector that avoids certain entropy violating solutions. Sargsyan et al.  have proposed a two-step method for uncertainty quantification of systems with discontinuous response. In the first step, a limited number of simulation runs are used to infer the shock location using the Bayesian inference. In the second step, localized gPC basis is defined in the region of analyticity, while intrusive Galerkin projection is used to propagate the uncertainty. In an alternate non gPC setting, Chantrasmi et al.  have proposed a Pade–Legendre interpolant-based approach, where simulation runs are obtained at predefined Gauss–Legendre quadrature nodes and resultant discontinuous system response is reconstructed using the Pade interpolation. Despite this recent progress, the resolution of Gibbs phenomenon for stochastic simulation of systems with discontinuous solution using SSP methods has not been fully resolved and is an area of current research interest . In addition, in the methods based on the gPC setting, the computational cost grows quickly as the stochastic dimension increase.
The rest of the paper is organized as follows: In Sec. 2, mathematical formulation of the method is provided. Proposed postprocessing method for Resolution of the Gibbs phenomenon using dynamic bi-orthogonality-based method is discussed in Sec. 3. Section 4 provides numerical results to demonstrate efficacy of the proposed method. Finally, paper is summarized and concluded in Sec. 5.
Dynamically Bi-Orthogonal Field Equations (DBFE)
Definition of the following operators will be used in the remainder of the paper.
where N is the number of eigenfunctions retained in the expansion (10). Note that the DO condition constrains the dynamic evolution of the eigenfunction basis such that the orthonormality of the eigenfield is retained at all time steps .
Using the DO condition and bi-orthogonal expansion of , SPDE (1) can be reformulated into a set of N + 1 PDEs and N × P ordinary differential equations (ODEs) as follows.
where Σ is the covariance matrix with ijth element .
Initial and Boundary Conditions for DBFE.
where Σ is a covariance matrix.
Resolution of the Gibbs Phenomena Using Dynamically Bi-Orthogonal Field Equations
Key challenges involved in the numerical implementation of the DBFE for SPDEs with discontinuous solutions are: (a) numerical evaluation of the mean and the eigenfunctions (using Eqs. (17) and (18)) results in the development of spurious oscillations as discontinuities evolve in the solution; (b) the resultant bi-orthogonal expansion (10) exhibits the Gibbs phenomenon characterized by the oscillations near the discontinuity location and the slow convergence away from the discontinuity location. This paper proposes a two-step approach for the application of DBFE to the stochastic systems with discontinuous solution that exploits the orthogonal decomposition of the spatial field in the DBFE approach. In the first step, extension of the existing schemes is proposed to derive the nonoscillatory numerical scheme for DBFE. In the second step, a Gegenbauer reprojection-based method is proposed to mitigate the effects of the Gibbs phenomenon. In this section, the proposed approach is described in detail.
Numerical Scheme for Nonoscillatory Solution.
Due to discontinuities, numerical solution of Eqs. (17)–(19) is expected to contain spurious oscillations with reduced accuracy in and . Total variation diminishing (TVD) finite difference schemes are widely used in the literature to obtain nonoscillatory accurate weak solutions to hyperbolic partial differential equations . This paper demonstrates extension of a TVD scheme for numerical solution of Eqs. (17)–(19); however, note that other numerical schemes, like essentially nonoscillatory schemes, Reimann solvers, Godunov approach, etc. (see Ref.  for details), can similarly be extended for DBFE.
where and .
The numerical flux (46) is used in Eq. (40) to obtain the first-order central difference scheme for numerical solution of Eqs. (17)–(19). Fourth-order Runge–Kutta method is used for time integration in Eqs. (17) and (18), while Euler method is used to solve ODE (19). Gauss–Legendre quadrature is used to calculate inner products in spatial dimensions, while Monte Carlo integration is used for inner products in stochastic dimensions.
Postprocessing Using Gegenbauer Polynomial.
Let the solution of the SPDE in Eq. (1) is obtained using the DBFE method till time t = T such that the discontinuities are developed. For a given ω, let denote an interval of analyticity for , where a random variable can be defined such that . Gibbs phenomenon can be resolved in the interval by reprojecting the partial sum (10) on a suitable Gibbs complementary polynomial basis . In the present paper, the partial sum is reprojected on a Gegenbauer polynomial, which is defined for a as follows.
The partial sum (56) approximates with exponential accuracy in the subinterval . Note that determination of domain requires location of discontinuity. An edge detection method [37,38] can be used to locate points of discontinuities.
Note that the proposed method is implemented in two steps. The first step involves time evolution of the mean, eigenvalues, and eigenfunctions of the DBFE expansion. The proposed extension of the TVD scheme ensures unphysical conditions are not developed during this time evolution. However, the presence of the discontinuity in the solution results in poor quantification of uncertainty due to Gibbs phenomenon. Reprojection on the Gegenbauer polynomials recovers accuracy of this uncertainty quantification. The numerical implementation of this proposed approach is summarized in the Algorithm 1.
Algorithm 1 Numerical Implementation
1: procedure Time Evolution of DBFE Expansion
2: Obtain DBFE expansion of initial and Boundary condition.
3: Use numerical flux given by Eq. (46) with the Runge–Kutta method to obtain time evolution of mean, eigenfunctions and gPC expansion coefficients.
4: At an arbitrary time instance T, use the DBFE expansion to obtain necessary statistics of the solution.
5: procedure Reprojection on Gegenbauer Polynomials
6: At an arbitrary time instance T, use expansion (56) to reproject the solution on Gegenbauer polynomials.
7: Obtain the DBFE expansion from the reprojected solution and restart the time evolution.
Stochastic Burgers Equation.
Efficacy of the proposed method is investigated for solution of a stochastic one-dimensional Burgers equation with uncertain input conditions. Note that previous work has utilized deterministic Burgers equation extensively , while there have been recent investigation on solution of stochastic Burgers equation [30,40]. Coupled with the ease of implementation, stochastic Burgers equation provides an appropriate context to investigate the efficacy of the proposed method for mitigation of the Gibbs phenomenon.
Figure 1 shows deterministic solution of the Burgers equation at different time steps for mean initial condition with ub = 0 and . For the deterministic numerical solution, explicit fourth-order Runge–Kutta scheme is used for time evolution while the first-order central Kurganov and Tadmor scheme  (see Eq. (44)) is used in spatial dimension. One-dimensional grid is used with while time step of is used. With time, left state moves toward increasing x, while at time , when the characteristic lines first cross, solution has infinite slope at x = 0, while the solution becomes discontinuous at x = 0 for . Final solution is shown at time .
Solution of Stochastic Burgers Equation Using DBFE.
The differential operator (62) is used in Eq. (46) to derive the DBFE governing equations (17)–(19). The field equations are spatially discretized with , while the fourth-order explicit Runge–Kutta method with time step is used for the time integration.
Initial condition for the mean is specified using Eq. (58), while the initial condition for is given by eigenfunctions of the covariance function (59). The initial conditions for expansion coefficients are specified as per Eq. (34). Note that the eigenvalues for the covariance function (59) decrease rapidly with only first three eigenmodes dominant ; thus, N = 3 suffices for the accurate approximation of the uncertain initial condition. However, to account for a highly nonlinear nature of the Burgers equation, higher value of N may be required. As can be inferred from Fig. 1, solution of stochastic Burgers equation is expected to evolve discontinuities. Effect of the Gibbs phenomenon for the resultant DBFE solution is mitigated by (a) obtaining nonoscillatory solution for mean and eigenfield and (b) using postprocessing on the resultant spectral expansion.
Figure 2 investigates efficacy of the proposed central KT scheme (see Sec. 3.1) to obtain nonoscillatory solution using the DBFE method for the eigenfield. The figure compares solution of the eigenfield obtained using three different numerical implementations of the DBFE method: (a) central difference scheme, where terms involving , M, and T are neglected in Eq. (46) (indicated as No Treatment in Fig. 2), (b) central KT scheme in mean that neglects terms involving and M in Eq. (46) (Mean Treatment in Fig. 2), and (c) central KT scheme (All Treatment in Fig. 2). With central difference scheme, oscillations are observed in the eigenfield. When central KT scheme is used only in mean field and central difference for eigenfield, oscillations are still present in eigenfield, though the amplitude is reduced, as can be observed for the first eigenfunction. When full central KT scheme is used, oscillations in eigenfunction are resolved.
To investigate the accuracy of the proposed method, the DBFE solution of the stochastic Burgers equation is compared with the Monte Carlo method using 10,000 samples. Figure 3(a) shows comparison of the 90% confidence bound for Monte Carlo and the DBFE method. The discrepancy in the confidence bound is due to the oscillations in the samples near the shock location, characterizing the Gibbs phenomenon for each sample. Amplitude of the oscillations is negligible for the samples near mean, which increases with the distance of the samples from the mean. Note that the distance of a sample from the mean is characterized by expansion coefficients and the associated gPC basis . For samples close to the mean, reconstruction of through the bi-orthogonal expansion (10) is dominated by , for which nonoscillatory solution is obtained using the numerical scheme (46). However, for samples away from the mean, effect of the Gibbs phenomenon becomes dominant through eigenfunctions . Figure 3(b) shows the L1-error of the DBFE method relative to the Monte Carlo method. error of the order of is observed away from the shock location; however, near the shock location, error of the order of is observed. Also, it is observed that the error does not reduce significantly with increasing N. From Figs. 3(a) and 3(b), it may be concluded that the L1-error is maximum at the shock location; and the amplitude of the maximum L1-error does not decrease by increasing the number of eigenfunctions, N. These observations are characteristics of the Gibbs phenomenon, indicating the need for appropriate postprocessing.
The Gegenbauer reprojection-based approach, outlined in Sec. 3.2, is used for postprocessing to mitigate the effect of the Gibbs phenomenon. The bi-orthogonal expansion of is reprojected on seventh-order Gegenbauer polynomials, while total seven coefficients are used in the expansion (56). Integrals involving the Gegenbauer polynomials are numerically evaluated using Gauss–Gegenbauer quadrature with 100 nodes, whereas the integrals involving the polynomial chaos basis are numerically solved using Monte Carlo integration. In the present paper, shock is assumed to be located at a point with highest slope where u = 0 line is crossed. However, the proposed method can be implemented using any edge detection method for shock localization.
Figure 3(c) shows comparison of 90% confidence bound of DBFE method after postprocessing with Monte Carlo simulation. Using postprocessing, agreement with Monte Carlo simulation is significantly improved. Figure 3(d) shows L1 error after postprocessing. Comparing with Fig. 3(b), it can be observed that the error is reduced at all locations by the order of . At mean shock location, L1 error is of the order of .
Figure 4 shows comparison of the mean and variance of the DBFE solution with the Monte Carlo simulation. The comparison is shown for both without and with postprocessing cases. For the higher moments, agreement with the Monte Carlo method improves significantly after the postprocessing. Mean of the solution obtained using the DBFE method without postprocessing matches well with the Monte Carlo method (shown in Fig. 4(a)). Note that the mean is not contaminated by the Gibbs phenomenon, resulting in the good match between DBFE method without postprocessing and the Monte Carlo method, while the subsequent postprocessing retains the accuracy of the mean. Figure 4(b) shows the variance as a function of spatial location x. The variance is low in the region away from the shock location, whereas significant increase in the variance is observed at mean shock location. Good agreement between the Monte Carlo and the DBFE method is observed in the region away from the shock location; however, nontrivial difference is observed without postprocessing near the mean shock location, with DBFE predicting higher variance than the Monte Carlo. The high variance in DBFE predictions results due to the oscillations near shock location characterizing the Gibbs phenomenon, which are mitigated using the postprocessing, providing the close agreement with the Monte Carlo method.
Computational Efficiency of the Proposed Method.
The computational cost of the numerical implementation of the proposed DBFE method is compared with the Monte Carlo and the gPC method (see Ref.  for details of the solution of stochastic Burgers equation using the gPC method). To investigate the computational efficiency of the proposed method, the Monte Carlo method, the gPC method, and the DBFE are implemented on a desktop computer with Intel Core i5 CPU. Total 10,000 samples are used for the Monte Carlo method, while third-order polynomial chaos basis is used for the DBFE and the gPC method. Time required for simulations is obtained using the fortran intrinsic routine cpu_time. The comparison of CPU time is shown in Fig. 5. Both the gPC and the DBFE methods provide considerable computational speed-up over the Monte Carlo method, with the computational time requirement increasing with the number of eigenfunctions used. Even for N = 7, speed-up of the order of 3 is obtained as compared to Monte Carlo method. The computational cost of the DBFE method is lower than the gPC method, with the ratio of the CPU time for gPC and the DBFE decreasing with increasing N. Although the postprocessing increases the computational cost of the DBFE method, the DBFE method still remains computationally efficient with the CPU time approaching that of the gPC with N. Note that the computational cost of the proposed method depends on the complexity of the differential operator ; however, the computational cost of the postprocessing remains constant irrespective of the complexity of the SPDE. The proposed DBFE method with the postprocessing provides the solution to SPDE with the accuracy comparable to the Monte Carlo method at a significantly lower computational cost.
Quasi-One-Dimensional Nozzle Flow.
where A is nozzle area, ρ is density, v is velocity, E is total energy, and p is the pressure. For notational convenience, define , and in the following.
Initial and Boundary Conditions.
The DBFE expansion for and is also initialized similarly.
Deterministic boundary conditions are specified at both the nozzle inlet and exit conditions. At the inlet boundary, density and temperature are specified, while the velocity is linearly interpolated from the adjacent grid points. Static pressure is specified at the exit boundary condition, whereas density and velocity are extrapolated from the adjacent grid points. Note that due to uncertain nozzle area, is uncertain despite the deterministic boundary conditions.
All the numerical results presented in this section are obtained using a first-order KT scheme in spatial dimension and fourth-order Runge–Kutta method for time integration, while spatiotemporal discretization is used with and . The results are shown after total time of . Figure 6(a) shows the mean nozzle area profile, while Fig. 6(b) shows corresponding deterministic solution. Existence of the shock discontinuity in the divergent region of the nozzle can be observed, resulting in the subsonic flow at the nozzle exit.
Uncertainty in the nozzle area is specified using a Gaussian process defined in Eq. (63) with the variance . The uncertain nozzle area is represented using a Karhunen–Loeve expansion with the mean nozzle profile as shown in Fig. 6(a), while N = 3 eigenfunctions are used for the spectral expansion. Resultant initial conditions for the DBFE solution are obtained as per the discussion in Sec. 4.3. Figure 7 shows the eigenfunctions at time t = 0 (i.e., initial condition) and (i.e., final solution). First eigenfunction is changed significantly for all , and ; however, second and third eigenfunctions are close to the initial conditions for and . However, note that all the eigenfunctions have evolved significantly over time for . From the governing equations for quasi-one-dimensional nozzle flow, it can be noted that the uncertain nozzle area profile primarily affects the momentum equation; thus, significant evolution in all the eigenfunctions is observed for .
To investigate efficacy of the proposed method, solution of stochastic quasi-one-dimensional nozzle flow is obtained using the Monte Carlo method. Relevant statistics is estimated by collecting 10,000 samples. Figure 8 shows comparison of mean velocity and mean pressure obtained using the Monte Carlo and the proposed method. As can be observed from Fig. 8(a), mean velocity and pressure obtained using the DBFE match closely with the Monte Carlo method. The relative L1-error in mean is shown in Fig. 8(b). The relative L1-error is of the order of , with the maximum error near the mean shock location for both the velocity and pressure. Note that for the reason of brevity, results only for velocity and pressure are shown in this section, though other parameters show similar trend.
Figure 9 shows comparison of variance obtained using the Monte Carlo and the DBFE method. The variance predicted by the DBFE method near the mean shock location is higher than the variance predicted by the Monte Carlo method for both velocity and pressure. As already observed for the Burgers equation, the higher variance in the DBFE method results from the oscillations near discontinuities associated with the Gibbs phenomenon. Effect of the Gibbs phenomenon is mitigated using the proposed postprocessing method, as can be observed from Fig. 9 through the close match in the variance between the Monte Carlo and the DBFE method with postprocessing.
This paper has proposed a dynamic bi-orthogonality-based spectral projection method for uncertainty quantification for systems with discontinuous solutions. The proposed approach is implemented in two steps: in the first step, input uncertainty is propagated to the system response using DBFE method, while in the second step, effect of the Gibbs phenomenon is mitigated by reprojecting the mean and the eigenfield on the Gegenbauer polynomials. Efficacy of the proposed method has been investigated for solution of a one-dimensional stochastic Burgers equation and stochastic quasi-one-dimensional flow through a convergent-divergent nozzle. The numerical results presented in this paper have demonstrated the ability of the proposed postprocessing method to mitigate the effect of the Gibbs phenomenon as discontinuities develop in the solution. Note that the proposed method does not require a priori knowledge of the shock location; thus, a generic implementation for a variety of applications can be achieved by extending any numerical scheme to the DBFE, as demonstrated in this paper. The DBFE method is found to be computationally efficient than the gPC method; thus, the method becomes an attractive alternative to the gPC for solution of SPDEs.
This work was supported by AFOSR Grant (FA9550-12-1-0313).
Throughout this paper, the operator “=” is used interchangeably to represent “equal to” and “approximately equal to.” Note that all the instances of “=” in spectral representation defines the “approximately equal to” operator.