{
"metadata": {
"name": "",
"signature": "sha256:a7d635f37a3a118f65c5e2a051fd27b84dc7b41ff7627418beacef1a7d78fd1c"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Advection Equation and the REA algorithm\n",
"by Mauricio J. Del Razo S. 2014\n",
"\n",
"In this notebook we will develop and explain the REA (reconstruct, evolve and average algorithm) for the advection equation. Since any homogeneous hyperbolic linear system of equations can be decoupled into a system of advection equations, this finite volume method algorithm is a fundamental block in the numerical solution of hyperbolic problems. The ideas presented with this algorithm will remain to be relevant for non-linear systems."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Decoupling into a system of advection equations\n",
"\n",
"Consider a system of one dimensional linear hyperbolic equations, written in its first order form,\n",
"\n",
"$$\n",
"\\bar{q}_t + \\mathbf{A} \\bar{q}_x = 0. \n",
"$$\n",
"\n",
"We can decompose the matrix $\\mathbf{A} = \\mathbf{R \\Lambda R^{-1}}$, as the product of the matrix of eigenvectors $\\mathbf{R}$ with the diagonal matrix of eigenvalues $\\mathbf{\\Lambda}$ and the inverse of $\\mathbf{R}$. Substituiting into the previus equation and multipling by $\\mathbf{R^{-1}}$ on the left side, we obtain\n",
"\n",
"$$\n",
"\\mathbf{R^{-1}} \\bar{q}_t + \\mathbf{\\Lambda R^{-1}} \\bar{q}_x = 0. \n",
"$$\n",
"\n",
"We can do the change of variable $\\bar{w} = \\mathbf{R^{-1}}\\bar{q}$, and as $\\mathbf{R^{-1}}$ is constant in time and space, we obtain,\n",
"\n",
"$$\n",
"\\bar{w}_t + \\mathbf{\\Lambda} \\bar{w}_x = 0,\n",
"$$\n",
"\n",
"where $\\mathbf{\\Lambda}$ is the diagonal matrix of eigenvalues. This yields a system of uncoupled advection equations, where the propagation speeds correspond to the eigenvalues of $\\mathbf{\\Lambda}$ or $\\mathbf{A}$. This means that if we can solve an advection equation numerically, we can solve the original system of hyperbolic equations. \n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. The REA algorithm for the advection equation\n",
"\n",
"Consider now a simple one dimensional advection equation for some quantity $q$,\n",
"\n",
"$$\n",
"q_t + c q_x = 0, \\\\[3mm]\n",
"q(x,0) = q_0(x).\n",
"$$\n",
"\n",
"It has solution $q(x,t) = q_0(x-ct)$, like shown below,\n",
"\n",
"\n",
"\n",
"We would like to develop a numerical method that solves this equation and conserves the quantity $q$. We will first need to divide the space into grid cells $\\mathcal{C_i}=(x_{i-1/2},x_{i+1/2})$ separated by $\\Delta x$ and consider the cell averages of $q$ to be,\n",
"\n",
"$$\n",
"Q_i^n = \\frac{1}{\\Delta x} \\int_{x_{i-1/2}}^{x_{i+1/2}} q(x,t_n) dx\n",
"$$\n",
"\n",
"Once this is decided, we can begin the REA algorithm (for a detailed explanation see [LeVeque 2002](http://depts.washington.edu/clawpack/book.html) ):\n",
"\n",
"\n",
"
\n",
" Reconstruct a piecewise polynomial function $\\bar{q}^n(x,t_n)$ from the cell averages $Q_i^n$. The simplest possibility is a piecewsie constant $\\bar{q}^n(x,t_n) = Q_i^n$. More complicated approximations could be useful for better accuracy.\n",
"\n",
"
\n",
"\n",
"
\n",
" Evolve the advection equation a time $\\Delta t$ with the reconstruction as initial data to obtain $\\bar{q}^{n+1}(x,t_{n+1})$\n",
"\n",
"
\n",
"\n",
"
\n",
" Average $\\bar{q}^{n+1}(x,t_{n+1})$ over the cells to obtain new cell averages,\n",
"$$\n",
"Q_i^{n+1} = \\frac{1}{\\Delta x} \\int_{x_{i-1/2}}^{x_{i+1/2}} \\bar{q}^{n+1}(x,t_{n+1}) dx\n",
"$$\n",
"\n",
"\n",
"
\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3. A python implementation \n",
"\n",
"In this section, we will implement the REA algorithm for the advection equation. The code is well commented, so it can be easily followed by the reader."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Import some libraries\n",
"import numpy as np\n",
"from scipy.integrate import quad\n",
"\n",
"# Set advection speed (c in PDE)\n",
"c = 10.0\n",
"\n",
"# Set spatial domain \n",
"xmin = 0\n",
"xmax = 15\n",
"\n",
"# Define q0(x), initial function to be advected\n",
"def q0(x):\n",
" qq = np.exp(-(x-2.5)*(x-2.5))\n",
" return qq\n",
"\n",
"# Create a function for the REA algorithm (parameter: number of cells and final time)\n",
"def REA(ncells,tfinal):\n",
" \n",
" # Create grid edges, space step and init time\n",
" x = np.linspace(xmin,xmax,ncells+1) # grid edges\n",
" dx = x[1] - x[0]\n",
" time = 0\n",
" \n",
" # Define time step size obeying CFL condition\n",
" dt = dx/(2.0*c) \n",
" \n",
" # Initialize solution q(x,t): q\n",
" q = np.zeros(ncells)\n",
" \n",
" # Calculate cell averages at time 0 (using q0 function)\n",
" Qavg = np.zeros(ncells)\n",
" for i in range(ncells):\n",
" # Integrate q0 function in the grid intervals\n",
" integral = quad(q0, x[i], x[i+1]) \n",
" Qavg[i] = integral[0]/dx\n",
" \n",
" # Initialiaze REA algorithm \n",
" while (time <= tfinal):\n",
" # RECONSTRUCT STEP\n",
" q = 1.0*Qavg\n",
" \n",
" # EVOLVE and AVERAGE STEP (move by dt and average)\n",
" time = time + dt\n",
" for i in range(ncells):\n",
" if (c>=0):\n",
" Qavg[i] = (c*dt*q[i-1] + (dx-c*dt)*q[i])/dx\n",
" else:\n",
" Qavg[i] = (c*dt*q[i+1] + (dx-c*dt)*q[i])/dx\n",
" \n",
" # If final time close adjust dt to match desired final time\n",
" if (time > tfinal):\n",
" dt = tfinal - (time - dt)\n",
" for i in range(ncells):\n",
" Qavg[i] = (c*dt*q[i-1] + (dx-c*dt)*q[i])/dx\n",
" q = Qavg\n",
" \n",
" return (x,q)\n",
"\n",
"# Function to calculate exact solution for comparison\n",
"def q_exact(tfinal):\n",
" # Calculate dense x\n",
" xx = np.linspace(xmin,xmax,1000)\n",
" # Obtain exact solution\n",
" qexact = q0(xx-c*tfinal)\n",
" return (xx,qexact)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 5
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#Import plotting libraries\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"from pylab import *\n",
"from IPython.html.widgets import interact\n",
"from ipywidgets import StaticInteract, RangeWidget\n",
"\n",
"# Create plotting function for REA and exact solution\n",
"def plot_REA(ncells, tfinal):\n",
" # Create Figure to plot\n",
" fig = plt.figure()\n",
" ax = fig.add_subplot(111) \n",
" \n",
" # Call REA algorithm and exact solution\n",
" x, qREA = REA(ncells,tfinal)\n",
" xx, qexact = q_exact(tfinal)\n",
" \n",
" # Transform REA into piecewise function for pretty plotting\n",
" qREA_pp = 0*xx\n",
" for i in range(len(x)-1):\n",
" for j in range(len(xx)):\n",
" if (xx[j] >= x[i] and xx[j] < x[i+1]): \n",
" qREA_pp[j] = qREA[i]\n",
" \n",
" # Plot solutions\n",
" ax.plot(xx,qREA_pp,'-b',linewidth=2)\n",
" ax.plot(xx,qexact,'-g',linewidth=2)\n",
" ax.axis([min(xx),max(xx),0,1.2])\n",
" return fig\n",
"\n",
"# Create animated solution\n",
"StaticInteract(plot_REA, ncells=RangeWidget(60,760,100), tfinal=RangeWidget(0,0.8,0.05))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"\n",
" \n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
"
\n",
" \n",
" ncells: \n",
" \n",
"tfinal: \n",
"
\n",
" "
],
"metadata": {},
"output_type": "pyout",
"prompt_number": 8,
"text": [
""
]
}
],
"prompt_number": 8
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The exact solution is shown in green, while the one obtained by the REA algorithm is blue. The numerical diffusion is evident on lower grid resolutions and still present for higher grid resolutions; however, it can be improved using higher order mehtods by using a linear approximation in the reconstruction step instead of a constant piecewise averages approximation."
]
}
],
"metadata": {}
}
]
}
![](\"images/REA01.png\")
\n",
"![](\"images/REA02.png\")
\n",
"![](\"images/REA03.png\")
\n",
"![](\"images/advection.png\")
\n",
"\n",
"We would like to develop a numerical method that solves this equation and conserves the quantity $q$. We will first need to divide the space into grid cells $\\mathcal{C_i}=(x_{i-1/2},x_{i+1/2})$ separated by $\\Delta x$ and consider the cell averages of $q$ to be,\n",
"\n",
"$$\n",
"Q_i^n = \\frac{1}{\\Delta x} \\int_{x_{i-1/2}}^{x_{i+1/2}} q(x,t_n) dx\n",
"$$\n",
"\n",
"Once this is decided, we can begin the REA algorithm (for a detailed explanation see [LeVeque 2002](http://depts.washington.edu/clawpack/book.html) ):\n",
"\n",
"