This notebook contains material from cbe67701-uncertainty-quantification; content is available on Github.
4.2 Advection-Diffusion-Reaction (ADR) Example¶
Noah Wamble (nwamble@nd.edu) June 22nd, 2020
This example was adapted from code from:
McClarren, Ryan G (2018). Uncertainty Quantification and Predictive Computational Science: A Foundation for Physical Scientists and Engineers, Chapter 4: Local Sensitivity Analysis Based on Derivative Approximations, Springer, https://doi.org/10.1007/978-3-319-99525-0_4
4.2.1 What do we want to know?¶
- How would the QoIs vary for small, expected perturbations in the inputs if given knowledge of the QoI at a particular value? To accomplish this, we are interested in the derivative of the QoI at a nominal input value (95).
4.2.2 Why do we want to know that?¶
This information will allow us to see which parameters have a larger effect on the QoI.
Knowing the sensitivities allows the researcher to narrow the focus to a smaller number of parameters and then apply more rigorous techniques to only these parameters.
If the variabilities in these distributions are small, it may be possible to quantify the uncertainty in the QoI to these parameters using only local information (95).
4.2.3 4.3.1 Simple ADR Example¶
We will study the steady advection-diffusion-reaction equation in one-spatial dimension with a spatially constant, but uncertain, diffusion coefficient, a linear reaction term, and a prescribed uncertain source.
$$ \nu \frac{du}{dx} - \omega \frac{d^2u}{dx^2} + \kappa(x)u = S(x) $$With initial condition: $$ u(0)=u(10)=0 $$ Where: $$ \kappa(x) = \begin{cases} \kappa_l & x \: \epsilon \: (5,7.5)\\ \kappa_h & \text{otherwise}\\ \end{cases} $$ $$ $$ $$ S(x)= qx(10-x) $$
With $\mu_v =10$, $\mu_\omega= 20$, Var$(\mu_v)=0.0723493$, Var$(\omega)= 0.3195214$
$\mu_{\kappa_h} = 2$, Var$( \kappa_h ) = 0.002778142$ , and $\mu_{k_l} = 0.1$, Var$( \kappa_l )= 8.511570 \times 10^{-6}$, $\mu_q = 1,$ Var$(q)= 7.062353 \times 10^{-4}$
We also prescribe that the parameters ordered as $(\nu,\omega,\kappa_l,\kappa_u,q)$ have the following correlation matrix:
$$ R= \begin{pmatrix} 1.00 & 0.10 & -0.05 & 0.00 & 0.00\\ 0.10 & 1.00 & -0.40 & 0.30 & 0.50 \\ 0.00 & 0.30 & 0.20 & 1.00 & -0.10 \\ 0.00 & 0.50 & 0.00 & -0.10 & 1.00 \\ \end{pmatrix} $$The QoI is total reaction rate: $$ Q= \int^{10}_{0} \kappa(x_ u) u(x)dx $$
Import Libraries
import matplotlib.pyplot as plt
import numpy as np
import scipy.sparse as sparse
import scipy.sparse.linalg as linalg
import scipy.integrate as integrate
import math
Numerical Method to Solve ADR Equation
def ADRSource(Lx, Nx, Source, omega, v, kappa):
"""
This function solves the diffusion equation with a generalized source.
Inputs
Lx, Nx: Used to make the size of the dx term
Source: Source function
omega: diffusion coefficient
v: velocity field
kappa: reaction coefficient equation
Outputs:
Solution: u(x) values
Q: Total reaction rate
"""
#Initialize A matrix that will be solved
A = sparse.dia_matrix((Nx,Nx),dtype="complex")
#Initialize size of dx step
dx = Lx/Nx
i2dx2 = 1.0/(dx*dx)
#fill diagonal of A
A.setdiag(2*i2dx2*omega + np.sign(v)*v/dx + kappa)
#fill off diagonals of A
A.setdiag(-i2dx2*omega[1:Nx] +
0.5*(1-np.sign(v[1:Nx]))*v[1:Nx]/dx,1)
A.setdiag(-i2dx2*omega[0:(Nx-1)] -
0.5*(np.sign(v[0:(Nx-1)])+1)*v[0:(Nx-1)]/dx,-1)
#solving equation A x = Source
Solution = linalg.spsolve(A,Source)
#Use trapezoidal integrater to find integral
Q = integrate.trapz(Solution*kappa,dx=dx)
return Solution, Q
Define Parameters
#Define desired spatial map and step
Lx = 10
Nx = 2000
dx = Lx/Nx
#Define functions that will be used by numerical solver
#source function
Source_func = lambda x, q: q*x*(10-x)
#reaction constants
kappa_func = lambda x, kappal, kappah: kappah + (kappal-kappah)*(x>5)*(x<7.5)
#Initialize these to the correct size
v_func = lambda x,v: v*np.ones(x.size)
omega_func = lambda x,omega: omega*np.ones(x.size)
#Define array for plotting purposes
xs = np.linspace(dx/2,Lx-dx/2,Nx)
# evaluate source and kappa
source = Source_func(xs, 1)
kappa = kappa_func(xs, 0.1, 2)
#nominal values and given variances
omega_nom = 20
omega_var = 0.3195214
v_nom = 10
v_var = 0.0723493
kappal_nom = 0.1
kappal_var = 8.511570e-6
kappah_nom = 2
kappah_var = 0.002778142
q_nom = 1
q_var = 7.062353e-4
vs = v_func(xs, v_nom)
# print(vs)
#Call Numerical Function above to solve equation
sol,Q = ADRSource(Lx, Nx, source, omega_func(xs, omega_nom), vs, kappa)
# print(Q)
#Plot solution
plt.plot(xs,sol)
plt.xlabel('x')
plt.ylabel('u(x)')
plt.title('The solution u(x) evaluated at mean value of uncertain parameters ')
plt.show()
Sensitivity Analysis
Sensitivity using Forward Difference Derivative Approximation:
$$ \frac{\partial Q}{\partial x_i} \bigg|_{\overline{x}} \approx \frac{Q(\overline{x} + \delta_i\hat{e_i}) - Q(\overline{x})}{\delta_i} $$Scaled Sensitivity Coefficient:
$$ \overline{x_i}\frac{\delta Q}{\delta x_i}\bigg|_\overline{x} $$Sensitivity Index: $$ \sigma_i\frac{\delta Q}{\delta x_i}\bigg|_\overline{x} $$
#v sensitivity
#finite difference step amount
delta = 1e-6
vs_pert = v_func(xs, v_nom*(1+delta))
#Calculate second point for finite difference
sol,Q_new = ADRSource(Lx, Nx, source, omega_func(xs, omega_nom), vs_pert, kappa)
#finite difference Sensitivity Calculation
sens_v = (Q_new - Q)/(v_nom*delta)
#Scaled Sensitivity
scaled_v = sens_v*v_nom
#Sensitivity Index
index_v = math.sqrt(v_var)*sens_v
#Print Commands
print('######################### Sensitivity Analysis #################')
print("v sensitivity: ", sens_v)
print("v scaled sensitivity: ", scaled_v)
print("v Sensitivity index: ", index_v)
print('\n')
#omega sensitivity
#Calculate second point for finite difference
sol,Q_new = ADRSource(Lx, Nx, source, omega_func(xs, omega_nom*(1+delta)), vs, kappa)
#Sensitivity calculation using finite difference
sens_omega = (Q_new - Q)/(omega_nom*delta)
#Scaled sensitivity Calc
scaled_v = sens_omega*omega_nom
#Sensitivity Index Calc
index_v = math.sqrt(omega_var)*sens_omega
#Print Commands
print("omega sensitivity: ", sens_omega)
print("omega scaled sensitivity: ", scaled_v)
print("omega sensitivity index: ", index_v)
print('\n')
#kappal sensitivity
#Calculate second point value for finite difference
sol,Q_new = ADRSource(Lx, Nx, source, omega_func(xs, omega_nom), vs, kappa_func(xs, kappal_nom*(1+delta), 2))
#Sensitivity Calculation using finite difference
sens_kappal = (Q_new - Q)/(kappal_nom*delta)
#Scaled Sensitivity Calc
scaled_v = sens_kappal*kappal_nom
#Sensitivity Index Calc
index_v = math.sqrt(kappal_var)*sens_kappal
#Print Commands
print("Kappa_l sensitivity: ", sens_kappal)
print("Kappa_l scaled sensitivity: ", scaled_v)
print("Kappa_l sensitivity index: ", index_v)
print('\n')
#kappah sensitivity
#Calculate second point value for finite difference
sol,Q_new = ADRSource(Lx, Nx, source, omega_func(xs, omega_nom), vs, kappa_func(xs, kappal_nom, kappah_nom*(1+delta)))
#Sensitivity Calculation using finite difference
sens_kappah = (Q_new - Q)/(kappah_nom*delta)
#Scaled Sensitivity Calc
scaled_v = sens_kappah*kappah_nom
#Sensitivity Index Calc
index_v = math.sqrt(kappah_var)*sens_kappah
#Print Commands
print("Kappa_h sensitivity: ", sens_kappah)
print("Kappa_h scaled sensitivity: ", scaled_v)
print("Kappa_h sensitivity index: ", index_v)
print('\n')
#q sensitivity
#Calculate second point value for finite difference
sol,Q_new = ADRSource(Lx, Nx, Source_func(xs, q_nom*(1+delta)), omega_func(xs, omega_nom), vs,kappa_func(xs, kappal_nom, kappah_nom))
sens_q = (Q_new - Q)/(q_nom*delta)
#Scaled Sensitivity Calc
scaled_v = sens_q*q_nom
#Sensitivity Index Calc
index_v = math.sqrt(q_var)*sens_q
#print commands
print("q sensitivity: ", sens_q)
print("q scaled sensitivity: ", scaled_v)
print("q sensitivity index: ", index_v)
print('\n')