In [1]:
%matplotlib inline
# plots graphs within the notebook
%config InlineBackend.figure_format='svg' # not sure what this does, may be default images to svg format
from IPython.display import Image
from IPython.core.display import HTML
def header(text):
raw_html = '<h4>' + str(text) + '</h4>'
return raw_html
def box(text):
raw_html = '<div style="border:1px dotted black;padding:2em;">'+str(text)+'</div>'
return HTML(raw_html)
def nobox(text):
raw_html = '<p>'+str(text)+'</p>'
return HTML(raw_html)
def addContent(raw_html):
global htmlContent
htmlContent += raw_html
class PDF(object):
def __init__(self, pdf, size=(200,200)):
self.pdf = pdf
self.size = size
def _repr_html_(self):
return '<iframe src={0} width={1[0]} height={1[1]}></iframe>'.format(self.pdf, self.size)
def _repr_latex_(self):
return r'\includegraphics[width=1.0\textwidth]{{{0}}}'.format(self.pdf)
class ListTable(list):
""" Overridden list class which takes a 2-dimensional list of
the form [[1,2,3],[4,5,6]], and renders an HTML Table in
IPython Notebook. """
def _repr_html_(self):
html = ["<table>"]
for row in self:
html.append("<tr>")
for col in row:
html.append("<td>{0}</td>".format(col))
html.append("</tr>")
html.append("</table>")
return ''.join(html)
font = {'family' : 'serif',
'color' : 'black',
'weight' : 'normal',
'size' : 18,
}
A numerical simulation starts with the identification of a conceptual model, often a set of partial differential equations and boundary conditions, that best define the system of interest. The mathematical nature of the conceptual model and the available computing power dictates the choice of numerical methods. However any numerical simulation should go through rigorous steps of verification and validation before any result can be used. If possible, uncertainty quantification should also be performed. This notebook takes you through the main steps of verification and validation using the Burgers equation:
$$ \frac{\partial u}{\partial t}=u\frac{\partial u}{\partial t}+\nu\frac{\partial^2 u}{\partial x^2}\,. $$
This equation is solver over $x\in[0,2\pi]$ with periodic boundary condition: $$ u(x=0)=u(x=2\pi)\,. $$ With the initial boundary condition: \begin{eqnarray} u &=& -\frac{2\nu}{\phi}\frac{\partial \phi}{\partial x}+4\\ \phi & = & \exp\left(-\frac{x^2}{4\nu}\right)+\exp\left(-\frac{(x-2\pi)^2}{4\nu}\right) \end{eqnarray} \begin{eqnarray} u &=& -\frac{2\nu}{\phi}\frac{\partial \phi}{\partial x}+4\\ \phi & = & \exp\left(-\frac{(x-4t)^2}{4\nu(t+1)}\right)+\exp\left(-\frac{(x-4t-2\pi)^2}{4\nu(t+1)}\right) \end{eqnarray} Here is the code for the analytical solution using sympy.
In [33]:
import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
from sympy.utilities.lambdify import lambdify
x, nu, t = sp.symbols('x nu t')
#phi = sp.exp(-(x-4*t)**2/(4*nu))+sp.exp(-(x-4*t-2*np.pi)**2/(4*nu))
phi = sp.exp(-(x-4*t)**2/(4*nu*(t+1.)))+sp.exp(-(x-4*t-2*np.pi)**2/(4*nu*(t+1)))
dphidx = phi.diff(x)
u = -2*nu*dphidx/phi+4.
ufunc = lambdify((t, x, nu),u)
print('d phi/dx=',dphidx)
n = 200
L = 2.*np.pi
dx = L/n
nu = 0.1
nt = 200
dt = dx*nu
T = nt*dt
mesh = np.linspace(dx/2.,L-dx/2,n)
t = 0.
T = 0.1
uexact = np.array([ufunc(T,x,nu) for x in mesh])
plt.plot(mesh,uexact)
plt.xlim(0,L)
plt.show
print(T)
The process of determining that a model implementation accurately represents the developer’s conceptual description of the model and the solution to the model. (AIAA G-077-1998)
The process of determining the degree to which a model is an accurate representation of the real world from the perspective of the intended uses of the model (AIAA G-077-1998)
To solve the Burger equation, we need to choose a discretization method, here finite difference and then choose numerical schemes for to calculate the three derivatives. The computational domain is as follows:
In [3]:
PDF('figures/1D-grid-u-F.pdf',size=(600,200))
Out[3]:
Fig 1. Sketch of the computational domain. The velocity is stored at the center of the cells.
To solve the Burger equation, we consider the following numerical schemes:
$$ \left.u\frac{\partial u}{\partial x}\right\vert_{i}\approx\frac{u_{i+\frac{1}{2}}\overline{u}_{i+\frac{1}{2}}-u_{i-\frac{1}{2}}\overline{u}_{i-\frac{1}{2}}}{\Delta x}+{\cal O}(\Delta x^n) $$
where $\overline{u}_{i+\frac{1}{2}}$ is the interpolated velocity at the face $i+1/2$ of the cell according to:$$ \begin{split} u_{i+\frac{1}{2}}\overline{u}_{i+\frac{1}{2}}&\approx\max\left(0,u_{i+1/2}\right)\left[(1-g_1+g_2)u_{i}+g_1u_{i+1}-g_2u_{i-1}\right]\\ &-\min\left(0,u_{i+1/2}\right)\left[(1-g_1+g_2)u_{i+1}+g_1u_{i}-g_2u_{i+2}\right]+{\cal O}(\Delta x^m) \end{split} $$
and $$ u_{i+\frac{1}{2}}=\frac{u_{i+1}+u_{i}}{2}\,. $$$$ \frac{\partial^2 u}{\partial x^2}\approx\frac{u_{i-1}-2u_i+u_{i+1}}{\Delta x^2}+{\cal O}(\Delta x^2) $$
The process of diffusion is governed by viscous stresses for momentum or heat flux for temperature. It quantifies the amount of velocity or temperature leaving the system through diffusion. The viscous stress or heat flux are proportional to the gradient of velocity or temperature. For instance, the heat flux is $$ F_{i+\frac{1}{2}}=-k\frac{T_{i+1}-T_{i}}{\Delta x} $$ The temperature variation within cell $i$ is therefore (see Lecture 0): $$ \frac{\partial T}{\partial t}\Delta x\Delta y + \left(-k\frac{T_{i+1}-T_{i}}{\Delta x}\right)\Delta y+\left(k\frac{T_{i}-T_{i-1}}{\Delta x}\right)\Delta x=0\,, $$ which reduces to $$ \frac{\partial T}{\partial t}=k\frac{T_{i-1}-2T_i+T_{i+1}}{\Delta x^2}\,. $$ When $\Delta x\rightarrow 0$, you should recognize the heat equation: $$ \frac{\partial T}{\partial t}=k\frac{\partial^2 T}{\partial x^2}\,. $$\begin{eqnarray} u_i^{n+1/2} &=& u_i^n-\frac{\Delta t}{2}\left(RHS_i^n\right)\\ u_i^{n+1} &=& u_i^n-\Delta t\left(RHS_i^{n+1/2}\right) \end{eqnarray}
where $$ RHS_i^n = -\left.\left(u\frac{\partial u}{\partial x}\right)\right\vert_i^n+\nu\left.\frac{\partial^2 u}{\partial x^2}\right\vert_i^n $$
In [44]:
Nschemes_advection = 4
Scheme = np.array(['CS','US1','US2','US3'])
g_1 = np.array([1./2.,0.,0.,3./8.])
g_2 = np.array([0.,0.,1./2.,1./8.])
def interpolation4advection(a,N,num_scheme):
imask = np.where(Scheme == num_scheme)
g1 = g_1[imask]
g2 = g_2[imask]
f=np.zeros(N+1,dtype='float64')
zero_array = np.zeros(N+1,dtype='float64')
sign_v = np.zeros(N+1,dtype='float64')
sign_v[1:N] = np.sign(0.5*(a[0:N-1]+a[1:N]))
sign_v[0] = np.sign(0.5*(a[N-1]+a[0]))
sign_v[N] = sign_v[0]
f[2:N-1] = np.maximum(zero_array[2:N-1],sign_v[2:N-1])*((1.-g1+g2)*a[1:N-2]+g1*a[2:N-1]-g2*a[0:N-3])\
-np.minimum(zero_array[2:N-1],sign_v[2:N-1])*((1.-g1+g2)*a[2:N-1]+g1*a[1:N-2]-g2*a[3:N])
f[1] = np.maximum(zero_array[1],sign_v[1])*((1.-g1+g2)*a[0]+g1*a[1]-g2*a[N-1]) \
-np.minimum(zero_array[1],sign_v[1])*((1.-g1+g2)*a[1]+g1*a[0]-g2*a[2])
f[0] = np.maximum(zero_array[0],sign_v[0])*((1.-g1+g2)*a[N-1]+g1*a[0]-g2*a[N-2]) \
-np.minimum(zero_array[0],sign_v[0])*((1.-g1+g2)*a[0]+g1*a[N-1]-g2*a[1])
f[N] = f[0]
f[N-1] = np.maximum(zero_array[N-1],sign_v[N-1])*((1.-g1+g2)*a[N-2]+g1*a[N-1]-g2*a[N-3]) \
-np.minimum(zero_array[N-1],sign_v[N-1])*((1.-g1+g2)*a[N-1]+g1*a[N-2]-g2*a[0])
return f
table = ListTable()
table.append(['Scheme', '$g_1$', '$g_2$'])
for i in range(4):
table.append([Scheme[i],g_1[i], g_2[i]])
table
Out[44]:
In [39]:
def diffusion(a,N,dx,nu):
diff = np.zeros(N,dtype='float64')
diff[1:N-1] = nu*(a[0:N-2]-2.*a[1:N-1]+a[2:N])/dx**2
diff[0] = nu*(a[N-1]-2.*a[0]+a[1])/dx**2
diff[N-1] = nu*(a[N-2]-2.*a[N-1]+a[0])/dx**2
return diff
In [40]:
def divergence(f,N,dz):
div = np.zeros(N,dtype='float64')
div[0:N] = (f[1:N+1]-f[0:N])/dx
return div
In [64]:
Lx = 2.*np.pi
Nx = 400
dx = L/Nx
nu = 0.07
num_scheme = 'CS'
simulation_time = 0.5
dt = 0.001
x_u = np.linspace(dx/2., Lx-dx/2, Nx,dtype='float64')
x_flux = np.linspace(0., Lx, Nx,dtype='float64')
u = np.zeros(Nx,dtype='float64') #velocity at cell center
u_old = np.zeros(Nx,dtype='float64')
u_face = np.zeros(Nx+1,dtype='float64') #velocity at faces
flux = np.zeros(Nx+1,dtype='float64') #advection flux (located on faces)
RHS = np.zeros(Nx,dtype='float64') #right-hand-side terms (advection+diffusion)
#initialization:
uini = np.array([ufunc(0.,x,nu) for x in x_u],dtype='float64')
u = uini
dt = np.amin([dx/np.max(np.abs(u)),0.2*dx**2/nu,dt]) #voodoo magic
print('time step= %0.5f' %dt)
T = np.dtype('float64')
T = 0. # time
Nrk = 2
rk_coef = np.array([0.5,1.],dtype='float64')
#rk_coef = np.array([1.],dtype='float64')
while T < simulation_time:
#print(T)
uold = u
for irk in range(Nrk):
u_face = interpolation4advection(u,Nx,num_scheme)
RHS[0:Nx] = -u[0:Nx]*divergence(u_face,Nx,dx)
RHS += diffusion(u,Nx,dx,nu)
u = uold + rk_coef[irk]*dt*RHS
T += dt
#print(T)
#print(u)
uexact = np.array([ufunc(T,x,nu) for x in x_u],dtype='float64')
plt.plot(x_u,u,label='sim.')
plt.plot(x_u,uini,label='ini')
plt.plot(x_u,uexact,label='exact')
plt.xlabel('$x$', fontdict = font)
plt.ylabel('$u$', fontdict = font)
plt.xlim(0,L)
plt.legend(loc=3, bbox_to_anchor=[0, 1],
ncol=3, shadow=True, fancybox=True)
plt.show()
In [78]:
from IPython.display import clear_output
Lx = 2.*np.pi
Nx = 400
dx = L/Nx
nu = 0.1
num_scheme = 'US3'
simulation_time = 10.0
Tplot = 0.01
dt_fixed = 0.001
x_u = np.linspace(dx/2., Lx-dx/2, Nx,dtype='float64')
x_flux = np.linspace(0., Lx, Nx,dtype='float64')
u = np.zeros(Nx,dtype='float64') #velocity at cell center
u_old = np.zeros(Nx,dtype='float64')
u_face = np.zeros(Nx+1,dtype='float64') #velocity at faces
flux = np.zeros(Nx+1,dtype='float64') #advection flux (located on faces)
RHS = np.zeros(Nx,dtype='float64') #right-hand-side terms (advection+diffusion)
#initialization:
uini = np.ones(Nx,dtype='float64')
u = uini
uforce = np.array([ufunc(0.,x,nu) for x in x_u],dtype='float64')
dt = np.amin([dx/np.max(np.abs(u)),0.2*dx**2/nu,dt]) #voodoo magic
print('time step= %0.5f' %dt)
T = np.dtype('float64')
T = 0. # time
Nrk = 2
rk_coef = np.array([0.5,1.],dtype='float64')
#rk_coef = np.array([1.],dtype='float64')
Tp = 0
while T < simulation_time:
#print(T)
uold = u
dt = np.amin([dx/np.max(np.abs(u)),0.2*dx**2/nu,dt_fixed])
for irk in range(Nrk):
u_face = interpolation4advection(u,Nx,num_scheme)
RHS[0:Nx] = -u[0:Nx]*divergence(u_face,Nx,dx)
RHS += diffusion(u,Nx,dx,nu) + 0.1*uforce
u = uold + rk_coef[irk]*dt*RHS
T += dt
Tp += dt
if (Tp >= Tplot):
plt.plot(x_u,u,label='sim.')
plt.xlabel('$x$', fontdict = font)
plt.ylabel('$u$', fontdict = font)
plt.xlim(0,L)
plt.ylim(0,5.)
plt.legend(loc=3, bbox_to_anchor=[0, 1],
ncol=3, shadow=True, fancybox=True)
plt.show()
clear_output(wait=True)
Tp = 0.
#print(T)
#print(u)
#uexact = np.array([ufunc(T,x,nu) for x in x_u],dtype='float64')
#plt.plot(x_u,u,label='sim.')
#plt.xlabel('$x$', fontdict = font)
#plt.ylabel('$u$', fontdict = font)
#plt.xlim(0,L)
#plt.legend(loc=3, bbox_to_anchor=[0, 1],
# ncol=3, shadow=True, fancybox=True)
#plt.show()
In [79]:
plt.plot(x_u,u,label='sim.')
plt.xlabel('$x$', fontdict = font)
plt.ylabel('$u$', fontdict = font)
plt.xlim(0,L)
plt.legend(loc=3, bbox_to_anchor=[0, 1],
ncol=3, shadow=True, fancybox=True)
plt.show()
In [ ]: