Uncertainty quantification in computational fluid dynamics

In this tutorial we will learn

Polynomial chaos and Monte Carlo

This tutorial raises a very important question. Why bother using polynomials for estimating moments? What exactly is the advantage? Moreover, are we guaranteed that we will converge to the Monte Carlo solution? The answer is a resounding yes! Infact this is precisely what Dongbin Xiu and George Karniandakis showed in their seminal paper [1]. As always we begin with some definitions: Parameter, Poly and Basis.

from equadratures import *
import numpy as np

For our model problem, let’s consider Rosenbrock’s function

\[f(x_1, x_2) = (1 - x_1)^2 + 100(x_1 - x_2^2)^2,\]

where we will assume that \(x_1\) and \(x_2\) are two uncertainties. We will assume that both parameters are Gaussians with \(\mu=1\) and \(\sigma=2\). Our objective is to compute the mean and variance in the output. We start by defining our computational model

def rosenbrock_fun(x):
        return (1 - x[0])**2 + 100*(x[1] - x[0]**2)**2

Next, we set the number of evaluation points in each direction. Lets opt for 7 points along each direction—more than sufficient to approximate the function exactly.

mu = 1
sigma = 2
variance = sigma**2
x1 = Parameter(distribution="Gaussian", shape_parameter_A=mu, shape_parameter_B=variance, order=4)
x2 = Parameter(distribution="Gaussian", shape_parameter_A=mu, shape_parameter_B=variance, order=4)
parameters = [x1, x2]

Now, we can set the problem up, compute the coefficients, and then ask Effective Quadratures to output the mean and the variance.

parameters = [x1, x2]
basis = Basis('Tensor grid')
uqProblem = Poly(parameters, basis, method='numerical-integration')
uqProblem.computeCoefficients(rosenbrock_fun)
uqProblem.set_model(rosenbrock_fun)
mean, variance = uqProblem.get_mean_and_variance()
print(mean, variance)

>> 6804.000000000022 476659232.0000047

Now, we compare these results with Monte Carlo.

large_number = 1000000
s = sigma * np.random.randn(large_number,2) + mu
f = np.zeros((large_number,1))
for i in range(0, large_number):
    f[i,0] = rosenbrock_fun([s[i,0], s[i,1]])
print np.mean(f), np.var(f)

>> 6813.941920252046 483131338.5544447

The results are very close! In fact the polynomial approximation results are exact, because Rosenbrock’s function is a polynomial itself!

But what order should we use? This is a tough question to answer without any apriori knowledge of the function we wish to obtain statistical moments from. We defer this question to the later tutorials, but will explore the effect of the order on accuracy. The plots below show the convergence in mean and variance with different number of samples.

../_images/tutorial_4_fig_a.png

Figure. Comparative convergence of the mean.

../_images/tutorial_4_fig_b.png

CFD with two aleatory uncertainties

The full source code for this tutorial can be found here.

References

  • Xiu, D., Karniandakis, G. E., (2002). The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations. SIAM Journal on Scientific Computing, 24(2), Paper