Methods for Uncertainty Propagation and Quantiﬁcation in Numerical Models

Uncertainty Quantification (UQ) is becoming a greater concern as capabilities of numerical codes is steadily increasing. Indeed, because of ever more complex physical models and accurate computational methods used, the need for a characterization of uncertainties associated to incomplete knowledge of models’ data (such as boundary conditions, forcing terms and various models’ constants) is crucial to assess uncertainty in models output, predictability, respective impacts of different uncertainty sources and support subsequent decision making processes. We solve UQ problems in a probabilistic framework through intrusive stochastic spectral projections. It amounts to the resolution of a large set of coupled equations, called the spectral problem. The objective of the project is to proposed parallelized computational strategies for the efficient resolution of the spectral problems. 1 Stochastic spectral methods Consider a modelM with random data D(θ ∈ Θ) defined on a probability space (Θ, Σ, dμ). The model output, S(θ), is also a random quantity defined on the same probability space and with a dependence structure with regards to the data prescribed by the modelM: M(S(θ); D(θ)) = 0 for almost any θ ∈ Θ. 1.1 Stochastic basis Let ξ(θ) = {ξ1, . . . , ξN} ∈ Ξ a finite set of independent random variables, defined on (Θ, Σ, dμ), with known probability density function p. We denote L2(Ξ, p) the space of second order random functionals in ξ, S ∈ L2(Ξ, p) : E[S] = ∫

Expectation of U : Variance of U : Ψ k are orthogonal with : 1 Support of Ψ k limited to an element : Stochastic multi-element method [Deb et al., 2001], [Wang and Karniadakis, 2005] Fully decouple the approximation problem for = elements 2 Hierarchical orthogonal Ψ k : Stochastic Multiwavelet method [olm et al., 2004, 2006, 2009] Couple the approximation problem, well suited for adaptive strategy (MRA) Introduction Spectral UQ Solution methods Application to conduction Application to natural convection Conclusion Application to spectral UQ [Ghanem & Spanos, 1991] Input-data parametrization

Model
We assume that for a.e. ξ ∈ Ξ, the problem M(S, D(ξ)) = 0 1 is well-posed, 2 has a unique solution and that the random solution S(ξ) ∈ L 2 (Ξ, p ξ ) : Spectral UQ Solution methods Application to conduction Application to natural convection Conclusion Application to spectral UQ [Ghanem & Spanos, 1991] Input-data parametrization

Solution expansion
Spectral UQ Solution methods Application to conduction Application to natural convection Conclusion Application to spectral UQ [Ghanem & Spanos, 1991] Input-data parametrization

Solution expansion
Knowledge of the spectral coefficients s k fully determine the random solution.
Makes explicit the dependence between D(ξ) and S(ξ).
Need efficient procedure(s) to compute the s k . Best approximation is defined by minimizing a (weighted) sum of squares of residuals :

Advantages/issues
Convergence with number of regression points m Selection of the regression points and "regressors" Ψ k Error estimate Introduction Spectral UQ Solution methods Application to conduction Application to natural convection Conclusion Non-Intrusive methods Non intrusive spectral projection : NISP Exploit the orthogonality of the basis : Computation of (P + 1) N-dimensional integrals

MC LHS QMC
Quadrature points ξ (i) and weights w Set of P + 1 coupled problems.

Implicitly account for modes' coupling
Often inherit properties of the deterministic model

Requires adaptation of deterministic solvers
Treatment of non-linearities. Consider the linear steady heat equation in an isotropic two-dimensional domain Ω, with boundary ∂Ω.
x ∈ Ω → u(x) ∈ R is the temperature field satisfying : where ν > 0 is the thermal conductivity and f ∈ L 2 (Ω) is a source term. We consider homogeneous Dirichlet and Neumann conditions over the respective portions Γ d and Γ n of the domain boundary

Introduction
Spectral UQ Solution methods Application to conduction Application to natural convection Conclusion

Stochastic conduction problem
Let V be the functionals space on Ω such that : where H 1 (Ω) is the Sobolev space of square integrable functionals whose first order derivatives are also square integrable. The variational problem is : where N is the set of nodes which are not lying on Γ d and Φ i (x) are the shape functions associated to these nodes.

Introduction
Spectral UQ Solution methods Application to conduction Application to natural convection Conclusion

Stochastic conduction problem
The Galerkin formulation in V h is : The problem can be recast as a system of linear equations    where n = Card(N ).
[a] is a (sparse) symmetric positive definite matrix.

Introduction
Spectral UQ Solution methods Application to conduction Application to natural convection Conclusion

Stochastic conduction problem
Consider the case of random conductivity and source term, defined on an abstract probability space (Θ, Σ, P) : It shows that the semi-discrete solution consists in n = Card(N ) random variables U i (θ). They satisfy i∈N j∈N We assume ν and F parameterized with N independent r.v.
Examples of parameterization will be shown later. The space of second-order random functionals in ξ is spanned by the Polynomial Chaos basis : where the Ψ i 's form a set of orthogonal multidimensional polynomials in ξ : Provided that ν and F are second-order quantities, they have orthogonal representations : The stochastic expansions are truncated to a finite polynomial order No. Different orders of truncation may be considered for the conductivity, source and solution. For simplicity, we use the same truncation order No. It corresponds to a stochastic approximation space

Introduction
Spectral UQ Solution methods Application to conduction Application to natural convection Conclusion

Stochastic conduction problem
The Galerkin problem is obtained by inserting the expansions of ν, F , U h and test functions V ∈ V h ⊗ S P into the variational form of the semi discrete stochastic problem. This results in : It involves n × (P + 1) deterministic quantities Introduction Spectral UQ Solution methods Application to conduction Application to natural convection Conclusion Stochastic conduction problem Denote u k := (u 1,k . . . u n,k ) t ∈ R n the vector of nodal values of the k-th stochastic mode of the solution. With this notation, the Galerkin problem becomes : Find u 0 , . . . , u P such that for all k = 0, . . . , where the matrix A l has for coefficients A l i,j and the vector The matrix blocks are given by : The system [A]u = B is called the spectral or Galerkin problem. Existence and uniqueness of the solution : Properties of the Galerkin system have been the focus of many works. e.g. [Babuska,2002], [Babuska, 2005], [Frauenfelder, 2005], [Matthies, 2005] For Dirichlet boundary conditions, the Galerkin system for stochastic elliptic problems has a unique solution provided that the random conductivity field satisfies some probabilistic (sufficient) conditions.
For the deterministic discretization with P − 1 finite-elements, these probabilistic conditions reduce to The random conductivity β is assumed to be log-normal, with unit median value β = 1 and coefficient of variation C ≥ 1.
β is parametrized with a unique normalized Gaussian variable ξ 1 (θ) so N = 1, and the PC basis is made of the one-dimensional Hermite polynomials.
The PC coefficients β k have closed form expressions [Ghanem, 1999] : Convergence with the expansion order β being log-normal, so is its inverse, and the expansion of 1/β is consequently given by : The spectrum of the numerical solution should decay as |σ β | k /k!.
Consider the random conductivity field defined as : ν 1 and ν 2 are two independent log-normal random variables with respective medians ν 1 and ν 2 , and coefficients of variation C 1 and C 2 .
Two normalized Gaussian variables ξ 1 and ξ 2 are used to parametrize the conductivity field.

Introduction
Spectral UQ Solution methods Application to conduction Application to natural convection Conclusion for k such that α k 3 + α k 4 + α k 5 = No. For the spatial discretization, two finite-element meshes will be considered: a coarse mesh with 4,326 elements and a fine mesh with 9,140 elements. The two meshes are plotted in Fig. 5.13. These two meshes will be used to analyze the effect of the spatial discretization on the statistics of the solution.
Mean and standard deviation of the temperature field: Figure 5.14 shows the mean of the temperature field over the computational domain, computed using the fine finite-element mesh. The highest expected temperature is reported for a point on the boundary n , located along the domain mid-line and opposite to the hole containing the coolant fluid. This finding may have been anticipated from physical considerations, assuming that the heat flux n is large enough to heat the neighborhood of the first hole to a temperature greater than T e (otherwise, the maximum temperature would be observed on e ). The effect of the insulating casing is also seen from the large temperature gradient in the casing domain on the side of the first (hot) hole.

Heat Equation 133
Fig. 5.14 Mean temperature field for the problem of Sect. 5.1.5. The computations are performed using No = 4 (126 spectral modes) and the fine finite-element mesh

Heat Equation 133
Fig. 5.14 Mean temperature field for the problem of Sect. 5.1.5. The computations are performed using No = 4 (126 spectral modes) and the fine finite-element mesh  surely greater than T c . Although not negligible, these events have low predicted probabilities. It is also seen from Fig. 5.17 that W h e appears well converged with regard to the spatial discretization.
The pdfs of W h c are shown in Fig. 5.18. It is seen that for the two meshes the fluxes W h c are always negative, denoting that the heat always flows from the exchanger to the cooling fluid, a result consistent with the problem settings, for which T c < T e and W n ≥ 0 almost surely. However, in contrast with the results reported for the fluxes W h e in Fig. 5.17, significant differences are observed between the approximations on the two meshes. Specifically, compared to the coarse mesh predictions, the fine mesh results yield a distribution of the flux W h c which is shifted toward the lower values with a larger variability (broader distribution). The difference in the mean values of the fluxes is about 5 units. A deeper analysis reveals that these differences are essentially due to the slow convergence of the flux boundary condition W h n toward its exact value W n , as h → 0. In fact, for the finite-element method used here, one expects a convergence of the solution u h as O(h 2 ), but a slower convergence in O(h) or less for the fluxes. Alternative formulations, such as a mixed one, would be better suited to achieve a more accurate approximation of the fluxes with larger convergence rates. It is however stressed that the slow convergence is due to the deterministic discretization techniques, and is in no way related to the stochastic spectral decomposition of the random solution. 2 The latter still exhibits an exponential convergence rate with the expansion order No for fixed h.
Conservation errors are analyzed in Fig. 5.19, which shows the pdf of error term

Heat Equation 137
Fig. 5.20 Probability density functions of η h c and η h e of problem in Sect. 5.1.5 for the coarse (left) and fine (right) finite-element meshes. Solutions with No = 4 (126 spectral modes) are used the mesh is refined and a better agreement between the two estimates η h c and η h e for the finer mesh. This again denotes the decay of the conservation error as the mesh is refined. Interestingly enough, we see that the random efficiency exceeds unity with some small but non-vanishing probability. This is consistent with our previous finding for the random flux W e where it was found to be positive for a non-vanishing Pdfs of T max , flux at "cold" boundary and of efficiencies. Analysis of the variance (ANOVA) Any functional G(ξ 1 , . . . , ξ N ) ∈ L 2 (Ξ, P Ξ ) has unique orthogonal hierarchical decomposition of the form (the so-called Sobol or Hoeffding's decomposition) Orthogonality : Analysis of the variability of the flux W c (second line of Table 5.2) shows that it is essentially due to the variability in W n (T W n ≈ 0.92), while uncertainty in the conductivities are seen to induce a low variability, though the uncertainty in the conductivity of the casing (ν 1 ) has a larger impact than the uncertainty in the conductivity of the exchanger (ν 2 ): T ν 2 ≈ 2T ν 1 . Variabilities in the boundary temperatures have even weaker impact. Again, this finding could have been anticipated since most of the heat is effectively extracted through the boundary c , so that one could have reasonably expected the uncertainty on W n to have the largest impact on W c . However, assessment of the relative impacts of the temperatures and conductivities would have been difficult to anticipate.
More interesting is the analysis of the variability in the flux W e (third line of Table 5.2) which, we recall, corresponds to the flux leaking from the casing to the 130 5 Detailed Elementary Applications

Fig. 5.11
Computational domain for the heat exchanger problem of Sect. 5.1.5. n is a Neumann boundary with flux n ; c and e are Dirichlet boundaries with temperatures T c and T n , respectively. The conductivity of the inner domain is ν 1 while the conductivity of the casing is ν 2

BC and solution representations
Galerkin projection