Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Yves Dubief, 2015. NSF for support via NSF-CBET award #1258697. The code for the solution of the Burgers equation is a copy of http://engineering101.org/1d-burgers-equation/

Lecture 2: Temporal Stability



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,
        }

Overview of the problem

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 x}=\nu\frac{\partial^2 u}{\partial x^2}\,. $$

This equation is solved 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 [2]:

    
import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
from sympy.utilities.lambdify import lambdify

def wave(x_phi,N,xo):
    phi = np.zeros(N)
    phi = (1.0+np.cos(x_phi-xo))/2.
    xmask = np.where(np.abs(x_phi-xo) > np.pi)
    phi[xmask] = 0.
    return phi

Lx = 10.*np.pi
Nx = 200
dx = Lx/Nx
xini = 2.*np.pi

x_phi = np.linspace(dx/2.,Lx-dx/2.,Nx,dtype='float64')
phi =np.zeros(Nx,dtype='float64')
phi = wave(x_phi,Nx,xini)
plt.plot(x_phi,phi,lw=2)
plt.xlabel('$x$', fontdict = font)
plt.ylabel('$\phi$', fontdict = font)
plt.xlim(0,Lx)
plt.show()



In [3]:

    
from IPython.display import clear_output
Lx = 20.*np.pi
Nx = 200
dx = Lx/Nx
Simulation_time = 40.
Tplot = 1.0
dt =0.1
xini = 2.*np.pi
x_phi = np.linspace(dx/2.,Lx-dx/2.,Nx,dtype='float64')
RHS = np.zeros(Nx,dtype='float64')
phi = np.zeros(Nx,dtype='float64')
phi_exact = np.zeros(Nx,dtype='float64')
phi = wave(x_phi,Nx,xini)
u = 1.
def compute_rhs_2(u,phi,N):
    r = np.zeros(Nx,dtype='float64')
    r[1:Nx-1] = (phi[2:Nx] - phi[0:Nx-2])/(2.*dx)
    r[0] = 0.
    r[Nx-1] = 0.
    r *= u
    return r
def compute_rhs_1(u,phi,N):
    r = np.zeros(Nx,dtype='float64')
    r[1:Nx] = (phi[1:Nx] - phi[0:Nx-1])/dx
    r[0] = 0.
    r[Nx-1] = 0.
    r *= u
    return r
T = 0.
Tp = 0.
while T < Simulation_time:
    #print(T)
    RHS = compute_rhs_2(u,phi,Nx)
    phi += -dt*RHS
    T += dt
    Tp += dt
    if (Tp >= Tplot):
        plt.plot(x_phi,phi,lw=2,label='simu')
        xt = xini + T*u
        phi_exact = wave(x_phi,Nx,xt)
        plt.plot(x_phi,phi_exact,lw=2,label='exact')
        plt.xlabel('$x$', fontdict = font)
        plt.ylabel('$\phi$', fontdict = font)
        plt.xlim(0,Lx)
        plt.legend(loc=3, bbox_to_anchor=[0, 1],
           ncol=3, shadow=True, fancybox=True)
        plt.show()
        clear_output(wait=True)
        Tp = 0.

\begin{eqnarray} \phi_0^{n+1} &=& \phi_0^n -\frac{u\Delta t}{\Delta x}\left(\frac{\phi_{1}^{n+1}+\phi_0^{n+!}}{2}-\phi_{-\frac{1}{2}}\right)\\ \phi_i^{n+1} &=& \phi_i^n -\frac{u\Delta t}{\Delta x}\left(\frac{\phi_{i+1}^{n+1}+\phi_i^{n+1}}{2}-\frac{\phi_i^{n+1}+\phi_{i-1}^{n+1}}{2}\right)\\ \phi_{N_x-1}^{n+1} &=& \phi_{N_x-1}^n -\frac{u\Delta t}{\Delta x}\left(\phi_{N_x-1}^{n+1}-\phi_{N_x-2}^{n+1}\right) \end{eqnarray}

\begin{eqnarray} \left(1-\frac{u\Delta t}{2\Delta x}\right)\phi_0^{n+1}+\frac{u\Delta t}{2\Delta x}\phi_{1}^{n+1} &=& \phi_0^n+\frac{u\Delta t}{\Delta x}\phi_{-\frac{1}{2}}\\ -\frac{u\Delta t}{2\Delta x}\phi_{i-1}^{n+1}+\phi_i^{n+1}+\frac{u\Delta t}{2\Delta x}\phi_{i+1}^{n+1} &=& \phi_i^n\\ -\frac{u\Delta t}{\Delta x}\phi_{N_x-2}^{n+1}+\left(1+\frac{u\Delta t}{\Delta x}\right)\phi_{N_x-1}^{n+1}+ &=& \phi_{N_x-1}^n \end{eqnarray}



In [36]:

    
def build_matrix(nx,dt,dx,u):
    A = np.zeros((nx,nx),dtype='float64')
    for i in range (1,nx-1):
            im = i - 1
            ip = i + 1
            A[i, i] = 1. 
            A[i, im] = -u*dt/(2.*dx)
            A[i, ip] = u*dt/(2.*dx)
            A[0, 0] = 1.-u*dt/(2.*dx)
            A[0, 1] = u*dt/(2.*dx)
            A[nx-1, nx-2] = -u*dt/dx
            A[nx-1, nx-1] = 1+u*dt/dx
    return A
            
from IPython.display import clear_output
Lx = 20.*np.pi
Nx = 200
dx = Lx/Nx
Simulation_time = 40.
Tplot = 1.0
dt =0.01
xini = 2.*np.pi
phi_bc_0 = 0.
x_phi = np.linspace(dx/2.,Lx-dx/2.,Nx,dtype='float64')
RHS = np.zeros(Nx,dtype='float64')
phi = np.zeros(Nx,dtype='float64')
phi_old = np.zeros(Nx,dtype='float64')
phi_exact = np.zeros(Nx,dtype='float64')
A = np.zeros((Nx,Nx),dtype='float64')
b = np.zeros(Nx,dtype='float64')
phi = wave(x_phi,Nx,xini)
u = 1.
A = build_matrix(Nx,dt,dx,u)
T = 0.
Tp = 0.
while T < Simulation_time:
    #print(T)
    b = phi
    b[0] += u*dt/dx*phi_bc_0
    phi = np.linalg.solve(A, b)
    T += dt
    Tp += dt
    if (Tp >= Tplot):
        plt.plot(x_phi,phi,lw=2,label='simu')
        xt = xini + T*u
        phi_exact = wave(x_phi,Nx,xt)
        plt.plot(x_phi,phi_exact,lw=2,label='exact')
        plt.xlabel('$x$', fontdict = font)
        plt.ylabel('$\phi$', fontdict = font)
        plt.xlim(0,Lx)
        plt.legend(loc=3, bbox_to_anchor=[0, 1],
           ncol=3, shadow=True, fancybox=True)
        plt.show()
        clear_output(wait=True)
        Tp = 0.



In [37]:

    
def build_matrix(nx,dt,dx,u):
    A = np.zeros((nx,nx),dtype='float64')
    for i in range (1,nx-1):
            im = i - 1
            ip = i + 1
            A[i, i] = 1. 
            A[i, im] = -u*dt/(2.*dx)
            A[i, ip] = u*dt/(2.*dx)
            A[0, 0] = 1.-u*dt/(2.*dx)
            A[0, 1] = u*dt/(2.*dx)
            A[nx-1, nx-2] = -u*dt/dx
            A[nx-1, nx-1] = 1+u*dt/dx
    return A
            
from IPython.display import clear_output
Lx = 20.*np.pi
Nx = 200
dx = Lx/Nx
Simulation_time = 40.
Tplot = 1.0
dt =0.01
xini = 2.*np.pi
phi_bc_0 = 0.
x_phi = np.linspace(dx/2.,Lx-dx/2.,Nx,dtype='float64')
RHS = np.zeros(Nx,dtype='float64')
phi = np.zeros(Nx,dtype='float64')
phi_old = np.zeros(Nx,dtype='float64')
phi_exact = np.zeros(Nx,dtype='float64')
Nrk = 2
rk_coef = np.array([1./2.,1.],dtype = 'float64')
A_0 = np.zeros((Nx,Nx),dtype='float64')
A_1 = np.zeros((Nx,Nx),dtype='float64')
b = np.zeros(Nx,dtype='float64')
phi = wave(x_phi,Nx,xini)
u = 1.
A_0 = build_matrix(Nx,rk_coef[0]*dt,dx,u)
A_1 = build_matrix(Nx,rk_coef[1]*dt,dx,u)
T = 0.
Tp = 0.
while T < Simulation_time:
    #print(T)
    phi_old = phi
    for irk in range(Nrk):
        b = phi_old
        b[0] += rk_coef[irk]*u*dt/dx*phi_bc_0
        if (irk == 0):
            phi = np.linalg.solve(A_0, b)
        else:
            phi = np.linalg.solve(A_1, b)
    T += dt
    Tp += dt
    if (Tp >= Tplot):
        plt.plot(x_phi,phi,lw=2,label='simu')
        xt = xini + T*u
        phi_exact = wave(x_phi,Nx,xt)
        plt.plot(x_phi,phi_exact,lw=2,label='exact')
        plt.xlabel('$x$', fontdict = font)
        plt.ylabel('$\phi$', fontdict = font)
        plt.xlim(0,Lx)
        plt.legend(loc=3, bbox_to_anchor=[0, 1],
           ncol=3, shadow=True, fancybox=True)
        plt.show()
        clear_output(wait=True)
        Tp = 0.

\begin{eqnarray} \phi_0^{n+1} &=& \phi_0^n -\frac{u\Delta t}{\Delta x}\left(\frac{\phi_{1}^{n+1}+\phi_0^{n+1}}{2}-\phi_{-\frac{1}{2}}^{n+1}\right)\\ \phi_i^{n+1} &=& \phi_i^n -\frac{u\Delta t}{\Delta x}\left(\frac{\phi_{i+1}^{n+1}+\phi_i^{n+1}}{2}-\frac{\phi_i^{n+1}+\phi_{i-1}^{n+1}}{2}\right)\\ \phi_{N_x-1}^{n+1} &=& \phi_{N_x-1}^n -\frac{u\Delta t}{\Delta x}\left(\phi_{N_x-1}^{n+1}-\phi_{N_x-2}^{n+1}\right) \end{eqnarray}



In [50]:

    
def build_matrix(nx,dt,dx,u):
    A = np.zeros((nx,nx),dtype='float64')
    for i in range (1,nx-1):
            im = i - 1
            ip = i + 1
            A[i, i] = 1. 
            A[i, im] = -u*dt/(2.*dx)
            A[i, ip] = u*dt/(2.*dx)
            A[0, 0] = 1.+u*dt/(2.*dx)
            A[0, 1] = u*dt/(2.*dx)
            A[nx-1, nx-2] = -u*dt/dx
            A[nx-1, nx-1] = 1+u*dt/dx
    return A
x,u,t = sp.symbols('x u t')
wave_exact = sp.sin(x-u*t)
phi_wave = lambdify((x,u,t),wave_exact)
from IPython.display import clear_output
Lx = 20.*np.pi
Nx = 512
dx = Lx/Nx
Simulation_time = 40.
Tplot = 1.0
dt =0.001
phi_bc_0 = 0.
x_phi = np.linspace(dx/2.,Lx-dx/2.,Nx,dtype='float64')
RHS = np.zeros(Nx,dtype='float64')
phi = np.zeros(Nx,dtype='float64')
phi_old = np.zeros(Nx,dtype='float64')
phi_exact = np.zeros(Nx,dtype='float64')
A = np.zeros((Nx,Nx),dtype='float64')
b = np.zeros(Nx,dtype='float64')
phi = np.array([phi_wave(x,u,0.) for x in x_phi], dtype = 'float64')
u = 1.
A = build_matrix(Nx,dt,dx,u)
T = 0.
Tp = 0.
while T < Simulation_time:
    #print(T)
    T += dt
    b = phi
    b[0] += u*dt/dx*phi_wave(0.,u,T)
    phi = np.linalg.solve(A, b)
    Tp += dt
    if (Tp >= Tplot):
        plt.plot(x_phi,phi,'r-',lw=2,label='simu')
        phi_exact = np.array([phi_wave(x,u,T) for x in x_phi], dtype = 'float64')
        plt.plot(x_phi,phi_exact,'b--',lw=2,label='exact')
        plt.xlabel('$x$', fontdict = font)
        plt.ylabel('$\phi$', fontdict = font)
        plt.xlim(0,Lx)
        plt.ylim(-1.1,1.1)
        plt.legend(loc=3, bbox_to_anchor=[0, 1],
           ncol=3, shadow=True, fancybox=True)
        plt.show()
        clear_output(wait=True)
        Tp = 0.

\begin{eqnarray} \phi_0^{n+1} &=& \phi_0^n -\frac{u\Delta t}{2\Delta x}\left(\left(\frac{\phi_{1}^{n+1}+\phi_0^{n+1}}{2}-\phi_{-\frac{1}{2}}^{n+1}\right)+\left(\frac{\phi_{1}^{n}+\phi_0^{n}}{2}-\phi_{-\frac{1}{2}}^{n}\right)\right)\\ \phi_i^{n+1} &=& \phi_i^n -\frac{u\Delta t}{2\Delta x}\left(\left(\frac{\phi_{i+1}^{n+1}+\phi_i^{n+1}}{2}-\frac{\phi_i^{n+1}+\phi_{i-1}^{n+1}}{2}\right)+\left(\frac{\phi_{i+1}^{n}+\phi_i^{n}}{2}-\frac{\phi_i^{n}+\phi_{i-1}^{n}}{2}\right)\right)\\ \phi_{N_x-1}^{n+1} &=& \phi_{N_x-1}^n -\frac{u\Delta t}{2\Delta x}\left(\left(\phi_{N_x-1}^{n+1}-\phi_{N_x-2}^{n+1}\right)+\left(\phi_{N_x-1}^{n}-\phi_{N_x-2}^{n}\right)\right) \end{eqnarray}

\begin{eqnarray} \left(1+\frac{u\Delta t}{4\Delta x}\right)\phi_0^{n+1}+\frac{u\Delta t}{4\Delta x}\phi_{1}^{n+1} &=& \phi_0^n+\frac{u\Delta t}{2\Delta x}\phi^{n+1}_{-\frac{1}{2}}-\frac{u\Delta t}{2\Delta x}\left(\frac{\phi_{1}^{n}+\phi_0^{n}}{2}-\phi_{-\frac{1}{2}}^{n}\right)\\ -\frac{u\Delta t}{4\Delta x}\phi_{i-1}^{n+1}+\phi_i^{n+1}+\frac{u\Delta t}{4\Delta x}\phi_{i+1}^{n+1} &=& \phi_i^n-\frac{u\Delta t}{2\Delta x}\left(\frac{\phi_{i+1}^{n}+\phi_i^{n}}{2}-\frac{\phi_i^{n}+\phi_{i-1}^{n}}{2}\right)\\ -\frac{u\Delta t}{2\Delta x}\phi_{N_x-2}^{n+1}+\left(1+\frac{u\Delta t}{2\Delta x}\right)\phi_{N_x-1}^{n+1} &=& \phi_{N_x-1}^n-\frac{u\Delta t}{2\Delta x}\left(\phi_{N_x-1}^{n}-\phi_{N_x-2}^{n}\right) \end{eqnarray}



In [10]:

    
def build_matrix(nx,dt,dx,u):
    A = np.zeros((nx,nx),dtype='float64')
    for i in range (1,nx-1):
            im = i - 1
            ip = i + 1
            A[i, i] = 1. 
            A[i, im] = -u*dt/(4.*dx)
            A[i, ip] = u*dt/(4.*dx)
            A[0, 0] = 1.+u*dt/(4.*dx)
            A[0, 1] = u*dt/(4.*dx)
            A[nx-1, nx-2] = -u*dt/(2.*dx)
            A[nx-1, nx-1] = 1+u*dt/(2.*dx)
    return A
x,u,t = sp.symbols('x u t')
wave_exact = sp.sin(x-u*t)
phi_wave = lambdify((x,u,t),wave_exact)
from IPython.display import clear_output
Lx = 20.*np.pi
Nx = 256
dx = Lx/Nx
Simulation_time = 100.
dt =0.01
Tplot = 1.0
Nplot = Simulation_time/Tplot
error = np.zeros(Nplot,dtype = 'float64')
energy = np.zeros(Nplot,dtype = 'float64')
energy_exact = np.zeros(Nplot,dtype = 'float64')
Time_plot = np.linspace(Tplot,Simulation_time,Nplot)

Nt = T/dt
phi_bc_0 = 0.
x_phi = np.linspace(dx/2.,Lx-dx/2.,Nx,dtype='float64')
RHS = np.zeros(Nx,dtype='float64')
phi = np.zeros(Nx,dtype='float64')
phi_old = np.zeros(Nx,dtype='float64')
phi_exact = np.zeros(Nx,dtype='float64')
A = np.zeros((Nx,Nx),dtype='float64')
b = np.zeros(Nx,dtype='float64')
phi = np.array([phi_wave(x,u,0.) for x in x_phi], dtype = 'float64')
u = 1.
A = build_matrix(Nx,dt,dx,u)
T = 0.
Tp = 0.
ip = 0
emask = np.where(np.abs(x_phi-(Lx-2.*np.pi)) <= np.pi )

while T < Simulation_time:
    #print(T)
    T += dt
    b = phi
    b[0] += u*dt/(2.*dx)*phi_wave(0.,u,T) - u*dt/(2.*dx)*((phi[0]+phi[1])/2.-phi_wave(0.,u,T-dt))
    b[1:Nx-1] -= u*dt/(4.*dx)*(phi[2:Nx]-phi[0:Nx-2])
    b[Nx-1] -= u*dt/(2.*dx)*(phi[Nx-1]-phi[Nx-2])
    phi = np.linalg.solve(A, b)
    Tp += dt
    if (Tp >= Tplot):
        plt.figure(figsize=(10, 8), dpi=160, facecolor='w', edgecolor='k')
        plt.subplot(2,2,1)
        plt.plot(x_phi,phi,'r-',lw=2,label='simu')
        phi_exact = np.array([phi_wave(x,u,T) for x in x_phi], dtype = 'float64')
        plt.plot(x_phi,phi_exact,'b--',lw=2,label='exact')
        plt.xlabel('$x$', fontdict = font)
        plt.ylabel('$\phi$', fontdict = font)
        plt.xlim(0,Lx)
        plt.ylim(-1.1,1.1)
        plt.legend(loc=3, bbox_to_anchor=[0, 1],
           ncol=3, shadow=True, fancybox=True)
        error[ip] =np.sqrt(np.sum(np.power(phi-phi_exact,2))/Nx)
        Time_plot[ip] = T
        energy[ip] = np.sum(np.power(phi[emask],2)*dx)/(2.*np.pi)
        energy_exact[ip] = np.sum(np.power(phi_exact[emask],2)*dx)/(2.*np.pi)
        plt.subplot(2,2,2)
        plt.plot(Time_plot[0:ip+1],error[0:ip+1],'r-',lw=2)
        plt.xlim(0,Simulation_time)
        plt.subplot(2,2,3)
        #plt.plot(x_phi,phi-phi_exact,'r-',lw=2,label='simu')
        plt.scatter(Time_plot[0:ip+1],energy[0:ip+1],marker='D',c='red',label='simu')
        plt.scatter(Time_plot[0:ip+1],energy_exact[0:ip+1],marker='o',c='blue',label='exact')
        plt.xlabel('$x$', fontdict = font)
        plt.ylabel('$\phi-\phi_{exact}$', fontdict = font)
        plt.xlim(0,Simulation_time)
        plt.legend(loc=3, bbox_to_anchor=[0, 1],
           ncol=3, shadow=True, fancybox=True)
        #plt.xlim(0,Lx)
        #plt.ylim(-1.1,1.1)
        plt.subplot(2,2,4)
        Nt = T/dt
        colors = x_phi
        plt.scatter(phi_exact,phi,c=colors)
        plt.xlabel('$\phi_{exact}$', fontdict = font)
        plt.ylabel('$\phi$', fontdict = font)
        plt.show()
        clear_output(wait=True)
        
        Tp = 0. 
        ip +=1



In [67]:

    
plt.plot(x_phi,phi-phi_exact,'r-',lw=2,label='simu')
plt.xlabel('$x$', fontdict = font)
plt.ylabel('$\phi-\phi_{exact}$', fontdict = font)
plt.xlim(0,Lx)
plt.ylim(-1.1,1.1)
#plt.legend(loc=3, bbox_to_anchor=[0, 1],
#           ncol=3, shadow=True, fancybox=True)
plt.show()



In [8]:

    
plt.figure(figsize=(10, 8), dpi=160, facecolor='w', edgecolor='k')
plt.subplot(2,2,1)
plt.plot(x_phi,phi,'r-',lw=2,label='simu')
phi_exact = np.array([phi_wave(x,u,T) for x in x_phi], dtype = 'float64')
plt.plot(x_phi,phi_exact,'b--',lw=2,label='exact')
plt.xlabel('$x$', fontdict = font)
plt.ylabel('$\phi$', fontdict = font)
plt.xlim(0,Lx)
plt.ylim(-1.1,1.1)
plt.legend(loc=3, bbox_to_anchor=[0, 1],
           ncol=3, shadow=True, fancybox=True)
plt.subplot(2,2,2)
plt.plot(Time_plot,error,'r-',lw=2)
plt.xlim(0,Simulation_time)
plt.xlabel('$t$', fontdict = font)
plt.ylabel('$\Vert\phi-\phi_{exact}\Vert_2$', fontdict = font)
plt.xlim(0,Simulation_time)
plt.subplot(2,2,3)
        #plt.plot(x_phi,phi-phi_exact,'r-',lw=2,label='simu')
plt.scatter(Time_plot,energy,marker='D',c='red',label='simu')
plt.scatter(Time_plot,energy_exact,marker='o',c='blue',label='exact')
plt.xlabel('$t$', fontdict = font)
plt.ylabel('energy', fontdict = font)
plt.xlim(0,Simulation_time)
#plt.legend(loc=3, bbox_to_anchor=[0, 1],
#           ncol=3, shadow=True, fancybox=True)
#plt.xlim(0,Lx)
#plt.ylim(-1.1,1.1)
plt.subplot(2,2,4)
colors = x_phi
plt.scatter(phi_exact,phi,c=colors)
plt.xlabel('$\phi_{exact}$', fontdict = font)
plt.ylabel('$\phi$', fontdict = font)
plt.savefig('figures/Lecture2-CK-CS2-1024-dt-01.pdf', bbox_inches='tight')
plt.show()



In [99]:

    
PDF('figures/Lecture2-CK-CS2-512-dt-001.pdf',size=(800,600))









    Out[99]:

\begin{eqnarray} \phi_0^{n+1} &=& \phi_0^n -\frac{u\Delta t}{2\Delta x}\left(\left(\phi_0^{n+1}-\phi_{-1}^{n+1}\right)+\left(\phi_0^{n}-\phi_{-1}^{n}\right)\right)\\ \phi_i^{n+1} &=& \phi_i^n -\frac{u\Delta t}{2\Delta x}\left(\left(\phi_i^{n+1}-\phi_{i-1}^{n+1}\right)+\left(\phi_i^{n}-\phi_{i-1}^{n}\right)\right)\\ \phi_{N_x-1}^{n+1} &=& \phi_{N_x-1}^n -\frac{u\Delta t}{2\Delta x}\left(\left(\phi_{N_x-1}^{n+1}-\phi_{N_x-2}^{n+1}\right)+\left(\phi_{N_x-1}^{n}-\phi_{N_x-2}^{n}\right)\right) \end{eqnarray}

\begin{eqnarray} \left(1+\frac{u\Delta t}{2\Delta x}\right)\phi_0^{n+1} &=& \phi_0^n+\frac{u\Delta t}{2\Delta x}\phi^{n+1}_{-1}-\frac{u\Delta t}{2\Delta x}\left(\phi_0^{n}-\phi_{-1}^{n}\right)\\ -\frac{u\Delta t}{2\Delta x}\phi_{i-1}^{n+1}+\left(1+\frac{u\Delta t}{2\Delta x}\right)\phi_i^{n+1} &=& \phi_i^n-\frac{u\Delta t}{2\Delta x}\left(\phi_i^{n}-\phi_{i-1}^{n}\right)\\ -\frac{u\Delta t}{2\Delta x}\phi_{N_x-2}^{n+1}+\left(1+\frac{u\Delta t}{2\Delta x}\right)\phi_{N_x-1}^{n+1} &=& \phi_{N_x-1}^n-\frac{u\Delta t}{2\Delta x}\left(\phi_{N_x-1}^{n}-\phi_{N_x-2}^{n}\right) \end{eqnarray}



In [13]:

    
def build_matrix_1(nx,dt,dx,u):
    A = np.zeros((nx,nx),dtype='float64')
    for i in range (1,nx-1):
            im = i - 1
            ip = i + 1
            A[i, i] = 1.+u*dt/(2.*dx) 
            A[i, im] = -u*dt/(2.*dx)
            A[i, ip] = 0.
            A[0, 0] = 1.+u*dt/(2.*dx)
            A[0, 1] = 0.
            A[nx-1, nx-2] = -u*dt/(2.*dx)
            A[nx-1, nx-1] = 1+u*dt/(2.*dx)
    return A
x,u,t = sp.symbols('x u t')
wave_exact = sp.sin(x-u*t)
phi_wave = lambdify((x,u,t),wave_exact)
from IPython.display import clear_output
Lx = 20.*np.pi
Nx = 512
dx = Lx/Nx
Simulation_time = 50.
dt =0.01
Tplot = 1.0
Nplot = Simulation_time/Tplot
error = np.zeros(Nplot,dtype = 'float64')
energy = np.zeros(Nplot,dtype = 'float64')
energy_exact = np.zeros(Nplot,dtype = 'float64')
Time_plot = np.linspace(Tplot,Simulation_time,Nplot)

Nt = T/dt
phi_bc_0 = 0.
x_phi = np.linspace(dx/2.,Lx-dx/2.,Nx,dtype='float64')
RHS = np.zeros(Nx,dtype='float64')
phi = np.zeros(Nx,dtype='float64')
phi_old = np.zeros(Nx,dtype='float64')
phi_exact = np.zeros(Nx,dtype='float64')
A = np.zeros((Nx,Nx),dtype='float64')
b = np.zeros(Nx,dtype='float64')
phi = np.array([phi_wave(x,u,0.) for x in x_phi], dtype = 'float64')
u = 1.
A = build_matrix_1(Nx,dt,dx,u)
T = 0.
Tp = 0.
ip = 0
emask = np.where(np.abs(x_phi-(Lx-2.*np.pi)) <= np.pi )

while T < Simulation_time:
    #print(T)
    T += dt
    b = phi
    b[0] += u*dt/(2.*dx)*phi_wave(-dx/2.,u,T) - u*dt/(2.*dx)*((phi[0]+phi[1])/2.-phi_wave(-dx/2.,u,T-dt))
    b[1:Nx-1] -= u*dt/(4.*dx)*(phi[2:Nx]-phi[0:Nx-2])
    b[Nx-1] -= u*dt/(2.*dx)*(phi[Nx-1]-phi[Nx-2])
    phi = np.linalg.solve(A, b)
    Tp += dt
    if (Tp >= Tplot):
        plt.figure(figsize=(10, 8), dpi=160, facecolor='w', edgecolor='k')
        plt.subplot(2,2,1)
        plt.plot(x_phi,phi,'r-',lw=2,label='simu')
        phi_exact = np.array([phi_wave(x,u,T) for x in x_phi], dtype = 'float64')
        plt.plot(x_phi,phi_exact,'b--',lw=2,label='exact')
        plt.xlabel('$x$', fontdict = font)
        plt.ylabel('$\phi$', fontdict = font)
        plt.xlim(0,Lx)
        plt.ylim(-1.1,1.1)
        plt.legend(loc=3, bbox_to_anchor=[0, 1],
           ncol=3, shadow=True, fancybox=True)
        error[ip] =np.sqrt(np.sum(np.power(phi-phi_exact,2))/Nx)
        Time_plot[ip] = T
        energy[ip] = np.sum(np.power(phi[emask],2)*dx)/(2.*np.pi)
        energy_exact[ip] = np.sum(np.power(phi_exact[emask],2)*dx)/(2.*np.pi)
        plt.subplot(2,2,2)
        plt.plot(Time_plot[0:ip+1],error[0:ip+1],'r-',lw=2)
        plt.xlim(0,Simulation_time)
        plt.subplot(2,2,3)
        #plt.plot(x_phi,phi-phi_exact,'r-',lw=2,label='simu')
        plt.scatter(Time_plot[0:ip+1],energy[0:ip+1],marker='D',c='red',label='simu')
        plt.scatter(Time_plot[0:ip+1],energy_exact[0:ip+1],marker='o',c='blue',label='exact')
        plt.xlabel('$x$', fontdict = font)
        plt.ylabel('energy', fontdict = font)
        plt.xlim(0,Simulation_time)
        #plt.legend(loc=3, bbox_to_anchor=[0, 1],
        #   ncol=3, shadow=True, fancybox=True)
        #plt.xlim(0,Lx)
        #plt.ylim(-1.1,1.1)
        plt.subplot(2,2,4)
        Nt = T/dt
        colors = x_phi
        plt.scatter(phi_exact,phi,c=colors)
        plt.xlabel('$\phi_{exact}$', fontdict = font)
        plt.ylabel('$\phi$', fontdict = font)
        plt.show()
        clear_output(wait=True)
        
        Tp = 0. 
        ip +=1



In [15]:

    
x,u,t = sp.symbols('x u t')
wave_exact = sp.sin(x-u*t)
phi_wave = lambdify((x,u,t),wave_exact)
from IPython.display import clear_output
Lx = 20.*np.pi
Nx = 128
dx = Lx/Nx
Simulation_time = 50.
dt =0.01
Tplot = 1.0
Nplot = Simulation_time/Tplot
error = np.zeros(Nplot,dtype = 'float64')
energy = np.zeros(Nplot,dtype = 'float64')
energy_exact = np.zeros(Nplot,dtype = 'float64')
Time_plot = np.linspace(Tplot,Simulation_time,Nplot)

Nt = T/dt
phi_bc_0 = 0.
x_phi = np.linspace(dx/2.,Lx-dx/2.,Nx,dtype='float64')
RHS = np.zeros(Nx,dtype='float64')
phi = np.zeros(Nx,dtype='float64')
phi_old = np.zeros(Nx,dtype='float64')
phi_exact = np.zeros(Nx,dtype='float64')
phi = np.array([phi_wave(x,u,0.) for x in x_phi], dtype = 'float64')
u = 1.
T = 0.
Tp = 0.
ip = 0
emask = np.where(np.abs(x_phi-(Lx-2.*np.pi)) <= np.pi )

while T < Simulation_time:
    #print(T)
    T += dt
    RHS[1:Nx] = -u*(phi[1:Nx]-phi[0:Nx-1])/dx
    RHS[0] = -u*(phi[0]-phi_wave(-dx/2.,u,T-dt))/dx
    phi = phi + dt*RHS
    Tp += dt
    if (Tp >= Tplot):
        plt.figure(figsize=(10, 8), dpi=160, facecolor='w', edgecolor='k')
        plt.subplot(2,2,1)
        plt.plot(x_phi,phi,'r-',lw=2,label='simu')
        phi_exact = np.array([phi_wave(x,u,T) for x in x_phi], dtype = 'float64')
        plt.plot(x_phi,phi_exact,'b--',lw=2,label='exact')
        plt.xlabel('$x$', fontdict = font)
        plt.ylabel('$\phi$', fontdict = font)
        plt.xlim(0,Lx)
        plt.ylim(-1.1,1.1)
        plt.legend(loc=3, bbox_to_anchor=[0, 1],
           ncol=3, shadow=True, fancybox=True)
        error[ip] =np.sqrt(np.sum(np.power(phi-phi_exact,2))/Nx)
        Time_plot[ip] = T
        energy[ip] = np.sum(np.power(phi[emask],2)*dx)/(2.*np.pi)
        energy_exact[ip] = np.sum(np.power(phi_exact[emask],2)*dx)/(2.*np.pi)
        plt.subplot(2,2,2)
        plt.plot(Time_plot[0:ip+1],error[0:ip+1],'r-',lw=2)
        plt.xlim(0,Simulation_time)
        plt.subplot(2,2,3)
        #plt.plot(x_phi,phi-phi_exact,'r-',lw=2,label='simu')
        plt.scatter(Time_plot[0:ip+1],energy[0:ip+1],marker='D',c='red',label='simu')
        plt.scatter(Time_plot[0:ip+1],energy_exact[0:ip+1],marker='o',c='blue',label='exact')
        plt.xlabel('$x$', fontdict = font)
        plt.ylabel('energy', fontdict = font)
        plt.xlim(0,Simulation_time)
        #plt.legend(loc=3, bbox_to_anchor=[0, 1],
        #   ncol=3, shadow=True, fancybox=True)
        #plt.xlim(0,Lx)
        #plt.ylim(-1.1,1.1)
        plt.subplot(2,2,4)
        Nt = T/dt
        colors = x_phi
        plt.scatter(phi_exact,phi,c=colors)
        plt.xlabel('$\phi_{exact}$', fontdict = font)
        plt.ylabel('$\phi$', fontdict = font)
        plt.show()
        clear_output(wait=True)
        
        Tp = 0. 
        ip +=1

Definition of Verification

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 code needs to be checked for programming errors.
Numerical errors and their impact on the conservation of mass and energy should be carefully assessed.
For steady problems, the solution needs to converge iteratively.
The solution should be independent of further resolution refinement of the computational grid.
All variables should have a consistent behavior.

Simulation of Burger equation

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:

Advection:
$$ \left.u\frac{\partial u}{\partial x}\right\vert_{i}\approx u_{i}\frac{\overline{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} \overline{u}_{i+\frac{1}{2}}&\approx\max\left(0,s_{i+1/2}\right)\left[(1-g_1+g_2)u_{i}+g_1u_{i+1}-g_2u_{i-1}\right]\\ &-\min\left(0,s_{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 $$ s_{i+\frac{1}{2}}=\begin{cases} 1 &\text{ if } \frac{u_{i+1}+u_{i}}{2}>0\\ -1 &\text{ if } \frac{u_{i+1}+u_{i}}{2}<0\\ 0 & \text{otherwise} \end{cases} \,. $$

upwind

Q1: In the case of CS and US3 (also called quick scheme), derive the truncation errors for the interpolation and the divergence of flux.

Diffusion:
$$ \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}\,. $$
Q2: Demonstrate that the numerical scheme used for the diffusion is second-order accurate in space.
Time variation of velocity: We will consider two schemes
1. The Euler scheme:
  \begin{eqnarray} u_i^{n+1} = u_i^n+\Delta t\left(RHS_i^n\right) \end{eqnarray}
2. The 2^nd order Runge Kutta scheme:
  \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 $$
Q3: Show that this scheme, the 2^nd Runge-Kutta, is second-order accurate in time.



In [4]:

    
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[4]:




Scheme $g_1$ $g_2$
CS 0.5 0.0
US1 0.0 0.0
US2 0.0 0.5
US3 0.375 0.125



In [5]:

    
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 [6]:

    
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 [44]:

    
from IPython.display import clear_output
Lx = 2.*np.pi
Nx = 500
dx = L/Nx
nu = 0.07 
num_scheme = 'US3'

simulation_time = 0.5
Tplot = 0.01 # period of plotting
dt_fixed = 0.001
Nrk = 2
if Nrk == 1:
    rk_coef = np.array([1.],dtype='float64') # Euler
elif Nrk == 2:
    rk_coef = np.array([0.5,1.],dtype='float64') # 2nd order Runge Kutta

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_fixed]) #voodoo magic
print('time step= %0.5f' %dt)
T = np.dtype('float64')
T = 0. # time
Tp = 0. # plotting time
print(T,simulation_time)
while T < simulation_time:
    uold = u
    if (np.amax(np.abs(u)) > 0.):
        dt = np.amin([dx/np.max(np.abs(u)),0.2*dx**2/nu,dt])
    else:
        dt = np.amin([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)
        u = uold + rk_coef[irk]*dt*RHS
    T += dt
    Tp += dt
    if (Tp >= Tplot):
        plt.plot(x_u,u,lw=2,label='simu')
        plt.plot(x_u,uini,lw=2,label='ini')
        uexact = np.array([ufunc(T,x,nu) for x in x_u],dtype='float64')
        plt.plot(x_u,uexact,lw=2,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()
        clear_output(wait=True)
        Tp = 0.



In [12]:

    
Lx = 2.*np.pi
simulation_time = 0.5
nu = 0.07 
dt_fixed = 0.0001
Nresolution = 5
resolution = np.array([50, 100, 250, 500, 1000],dtype='int32')
deltax = np.zeros(2)
deltax = Lx/resolution
print(deltax)

L2norm = np.zeros((Nresolution,Nschemes_advection,2))

for Nrk in np.array([1, 2]):
    print(Nrk)
    if Nrk == 1:
        rk_coef = np.array([1.],dtype='float64') # Euler
    elif Nrk == 2:
        rk_coef = np.array([0.5,1.],dtype='float64') # 2nd order Runge Kutta
    ires = -1
    for Nx in resolution:
        #print(Nx)
        ires += 1
        dx = L/Nx
        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)
        ischeme = -1
        for num_scheme in Scheme:
            #print(num_scheme)
            ischeme += 1
            #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_fixed]) #voodoo magic
            T = np.dtype('float64')
            T = 0. # time
            while T < simulation_time:
                uold = u
                if (np.amax(np.abs(u)) > 0.):
                    dt = np.amin([dx/np.max(np.abs(u)),0.2*dx**2/nu,dt])
                else:
                    dt = np.amin([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)
                    u = uold + rk_coef[irk]*dt*RHS
                T += dt
            uexact = np.array([ufunc(T,x,nu) for x in x_u],dtype='float64')
            #L2norm[ires,ischeme,Nrk] = np.linalg.norm(u-uexact)
            L2norm[ires,ischeme,Nrk-1] = np.sqrt(np.sum(np.power(u-uexact,2))/Nx)
print('Done')









    



[ 0.12566371  0.06283185  0.02513274  0.01256637  0.00628319]
1
2
Done



In [34]:

    
Nrk = 1
plt.loglog(deltax,L2norm[:,0,Nrk-1],'b--',lw=2,label='Euler-CS')
plt.loglog(deltax,L2norm[:,1,Nrk-1],'k--',lw=2,label='Euler-US1')
plt.loglog(deltax,L2norm[:,2,Nrk-1],'g--',lw=2,label='Euler-US2')
plt.loglog(deltax,L2norm[:,3,Nrk-1],'r--',lw=2,label='Euler-US3')
plt.xlabel('$\Delta x$', fontdict = font)
plt.ylabel('$\Vert u-u_{exact}\Vert_2$', fontdict = font)
#plt.xlim(0,L)
plt.legend(loc=3, bbox_to_anchor=[0, 1],
           ncol=3, shadow=True, fancybox=True)
plt.show()



In [25]:

    
Nrk = 2
plt.loglog(deltax,L2norm[:,0,Nrk-1],'b-',lw=2,label='RK2-CS')
plt.loglog(deltax,L2norm[:,1,Nrk-1],'k-',lw=2,label='RK2-US1')
plt.loglog(deltax,L2norm[:,2,Nrk-1],'g-',lw=2,label='RK2-US2')
plt.loglog(deltax,L2norm[:,3,Nrk-1],'r-',lw=2,label='RK2-US3')
plt.xlabel('$\Delta x$', fontdict = font)
plt.ylabel('$\Vert u-u_{exact}\Vert_2$', fontdict = font)
#plt.xlim(0,L)
plt.legend(loc=3, bbox_to_anchor=[0, 1],
           ncol=2, shadow=True, fancybox=True)
plt.show()



In [29]:

    
Nrk = 1
plt.loglog(deltax,L2norm[:,0,Nrk-1],'b--',lw=2,label='Euler-CS')
plt.loglog(deltax,L2norm[:,1,Nrk-1],'k--',lw=2,label='Euler-US1')
plt.loglog(deltax,L2norm[:,2,Nrk-1],'g--',lw=2,label='Euler-US2')
plt.loglog(deltax,L2norm[:,3,Nrk-1],'r--',lw=2,label='Euler-US3')
Nrk = 2
plt.loglog(deltax,L2norm[:,0,Nrk-1],'b-',lw=2,label='RK2-CS')
plt.loglog(deltax,L2norm[:,1,Nrk-1],'k-',lw=2,label='RK2-US1')
plt.loglog(deltax,L2norm[:,2,Nrk-1],'g-',lw=2,label='RK2-US2')
plt.loglog(deltax,L2norm[:,3,Nrk-1],'r-',lw=2,label='RK2-US3')
plt.xlabel('$\Delta x$', fontdict = font)
plt.ylabel('$\Vert u-u_{exact}\Vert_2$', fontdict = font)
#plt.xlim(0,L)
plt.legend(loc=3, bbox_to_anchor=[0, 1],
           ncol=2, shadow=True, fancybox=True)
plt.show()



In [15]:

    
plt.plot(x_u,u,label='simu',lw=2)
#plt.plot(x_u,uini,label='ini',lw=2)
uexact = np.array([ufunc(T,x,nu) for x in x_u],dtype='float64')
plt.plot(x_u,uexact,'r--',label='exact',lw=2)
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 [38]:

    
plt.plot(x_u,(u-uexact)/uexact,label='error',lw=2)
plt.xlabel('$x$', fontdict = font)
plt.ylabel('$\Vert u-u_{exact}\Vert$', fontdict = font)
plt.xlim(0,L)
plt.legend(loc=3, bbox_to_anchor=[0, 1],
           ncol=3, shadow=True, fancybox=True)
plt.show()



In [18]:

    
table = ListTable()
table.append(['Scheme', 'CS', 'US1', 'US2', 'US3'])
ires = 50
Nrk = 2
for ires in range(Nresolution):
    table.append([resolution[ires],L2norm[ires,0,Nrk-1], L2norm[ires,1,Nrk-1],L2norm[ires,2,Nrk-1],L2norm[ires,3,Nrk-1]])
table









    Out[18]:




Scheme CS US1 US2 US3
50 0.447754654442 0.757306270684 0.480387118707 0.518058693828
100 0.112206946975 0.632783607746 0.246014211548 0.203550200629
250 0.0146265974871 0.448379995725 0.0474965492007 0.0227969640314
500 0.00347246973913 0.301444126703 0.00787928329122 0.00329566458324
1000 0.000854762844227 0.181003916944 0.00135658185627 0.000503893630791



In [ ]:

Scheme	CS	US1	US2	US3
50	0.447754654442	0.757306270684	0.480387118707	0.518058693828
100	0.112206946975	0.632783607746	0.246014211548	0.203550200629
250	0.0146265974871	0.448379995725	0.0474965492007	0.0227969640314
500	0.00347246973913	0.301444126703	0.00787928329122	0.00329566458324
1000	0.000854762844227	0.181003916944	0.00135658185627	0.000503893630791

Lecture 2: Temporal Stability

Overview of the problem

Definition of Verification

Simulation of Burger equation

Q1: In the case of CS and US3 (also called quick scheme), derive the truncation errors for the interpolation and the divergence of flux.

Q2: Demonstrate that the numerical scheme used for the diffusion is second-order accurate in space.

Q3: Show that this scheme, the 2nd Runge-Kutta, is second-order accurate in time.

Q3: Show that this scheme, the 2^nd Runge-Kutta, is second-order accurate in time.