{
"metadata": {
"name": ""
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Evaluating Burgers Equation with different CFD Schemes\n",
"===\n",
"\n",
"We've already examined Burgers equation in both 1D and 2D ([Step 4](http://nbviewer.ipython.org/urls/github.com/barbagroup/CFDPython/blob/master/lessons/05_Step_4.ipynb) and [Step 8](http://nbviewer.ipython.org/urls/github.com/barbagroup/CFDPython/blob/master/lessons/10_Step_8.ipynb), respectively). Here we are going to restrict ourselves to the 1D Burgers equation and examine the role that different schemes play in discretizing non-linear 1st order hyperbolic equations. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider the 1D Burgers equation:\n",
"\n",
"$$\\frac{\\partial u}{\\partial t} = -u \\frac{\\partial u}{\\partial x}$$\n",
"\n",
"We want to represent this in conservative form so that we can better deal with potential shocks, which gives us:\n",
"\n",
"$$\\frac{\\partial u}{\\partial t} = - \\frac{\\partial }{\\partial x}(\\frac{u^2}{2})$$\n",
"\n",
"We can also write this as:\n",
"\n",
"$$\\frac{\\partial u}{\\partial t} = - \\frac{\\partial E}{\\partial x}$$\n",
"\n",
"if we take $E = \\frac{u^2}{2}$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Initial Conditions\n",
"-----\n",
"For each scheme, the initial conditions are \n",
"\n",
"$$u(x,0) = \\left\\{ \\begin{array}{cc}\n",
"1 & 0 \\leq x < 2 \\\\\n",
"0 & 2 \\leq x \\leq 4 \\\\ \\end{array} \\right.$$\n",
"\n",
"Assignment\n",
"----\n",
"First investigate the behavior of the Burgers function with $\\Delta t = \\Delta x = 1$, that is, with a Courant number of $\\frac{\\Delta t}{\\Delta x}=1$, then change the Courant number to $\\frac{\\Delta t}{\\Delta x}=0.5$. \n",
"\n",
"Also experiment with different mesh sizes to see what kind of behavior each scheme exhibits under different circumstances. \n",
"\n",
"Helper Code\n",
"---\n",
"Start by importing `numpy` and `matplotlib.pyplot`. The initial conditions are also defined below, along with a few helper functions to return values that are frequently required. "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%matplotlib inline\n",
"import numpy\n",
"from matplotlib import pyplot\n",
"\n",
"\n",
"#Basic initial condition parameters\n",
"#defining grid size, time steps, CFL condition, etc...\n",
"nx = 81\n",
"nt = 70\n",
"sigma = 1\n",
"dx = 4.0/nx\n",
"dt = sigma*dx\n",
"\n",
"#Define a quick function to set up the initial square wave condition\n",
"def u_ic():\n",
" u = numpy.ones(nx)\n",
" u[(nx-1)/2:]=0\n",
" return u\n",
"\n",
"#Define two lambda functions to help with common operations encountered\n",
"#with the Euler equations\n",
"utoE = lambda u: (u/2)**2\n",
"utoA = lambda u: u**2"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##Lax-Friedrichs\n",
"\n",
"Lax-Friedrichs is an explicit, 1st order scheme, using forward-difference in time and central-difference in space. However, notice that the $u^n_i$ value is calculated as the average of its adjacent cells. \n",
"\n",
"$$\\frac{u_i^{n+1}-\\frac{1}{2}(u^n_{i+1}+u^n_{i-1})}{\\Delta t} = -\\frac{E^n_{i+1}-E^n_{i-1}}{2 \\Delta x}$$\n",
"\n",
"Transposing this to solve for $u^{n+1}_i$ yields:\n",
"\n",
"$$u_i^{n+1} = \\frac{1}{2}(u^n_{i+1}+u^n_{i-1}) - \\frac{\\Delta t}{2 \\Delta x}(E^n_{i+1}-E^n_{i-1})$$\n",
"\n",
"This scheme is implemented in the function below. All the schemes in this notebook are wrapped in their own functions to help with displaying animations of the results. This is also good coding practice!\n",
"\n",
"In order to display our animations, we're going to hold the results of each timestep in `u` in a 2D array. So our array `un` will have `nt` rows and `nx` columns."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def laxfriedrichs(u, nt, dt, dx):\n",
" #initialize our results array with dimensions nt by nx\n",
" un = numpy.zeros((nt,len(u))) \n",
" #copy the initial u array into each row of our new array\n",
" un[:,:] = u.copy() \n",
" \n",
" '''\n",
" Now, for each timestep, we're going to calculate u^n+1, \n",
" then set the value of u equal to u^n+1 so we can calculate \n",
" the next iteration. For every timestep, the entire vector\n",
" u^n is saved in a single row of our results array un.\n",
" '''\n",
" for i in range(1,nt):\n",
" E = utoE(u)\n",
" un[i,1:-1] = .5*(u[2:]+u[:-2]) - dt/(2*dx)*(E[2:]-E[:-2])\n",
" un[i,0] = 1\n",
" u = un[i].copy()\n",
" \n",
" return un"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The above cell just defines the function which can execute the Lax-Friedrichs scheme, now it needs to be called. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##Lax Friedrichs Test 1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's try to display the results of our Lax Friedrichs scheme. First, we need to generate the results by calling our function `laxfriedrics`. "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"u = u_ic() #make sure that u is set to our expected initial conditions\n",
"un = laxfriedrichs(u,nt,dt,dx)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 3
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can ask NumPy to tell us the shape of our new array `un` and see if we have the dimensions we're expecting. We should have `nt` rows and `nx` columns. "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"un.shape"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 4,
"text": [
"(70, 81)"
]
}
],
"prompt_number": 4
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Looks good! Now in order to view the results in an animation, we'll import a few extra libraries. The `JSAnimation` library is not part of the standard IPython install (yet!) but you can click [here](http://dummy) for instructions on how to install it."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from matplotlib import animation\n",
"from JSAnimation.IPython_display import display_animation"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 5
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To view a nice interactive animation of our results, we'll need to define a few parameters.\n",
"\n",
"* First, we need a `figure`, which we'll also make a little larger than the default size for easier viewing. \n",
"* We also want to define an axis in the figure with some plotting limits. \n",
"* Then, we want to initialize a `line` object that will be a part of our plot, but we don't need to give it any data yet, so we'll just assign it empty lists for its `x` and `y` coordinates.\n",
"\n",
"Now we need a functions to \"draw\" each frame of our animation:\n",
"\n",
"* `animate` takes a single row of our results array `un` and then plots a line.\n",
"\n",
"Now to set up the `animation` object, we pass is a few values:\n",
"\n",
"* `fig` is the figure we've already defined and we're telling `FuncAnimation` to use that to draw in\n",
"* `animate` is the function we're using to actually plot our results\n",
"* `frames=un` tells `FuncAnimation` that it should use each row of our array `un` to generate sequential frames of the animation.\n",
"* `interval=50` sets the default amount of time in milliseconds between frames\n",
"\n",
"You'll notice a set of controls below the animation -- you can use the **+** and **-** buttons to speed up or slow down the animation."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"fig = pyplot.figure();\n",
"ax = pyplot.axes(xlim=(0,4),ylim=(-.5,2));\n",
"line, = ax.plot([],[],lw=2);\n",
"\n",
"def animate(data):\n",
" x = numpy.linspace(0,4,nx)\n",
" y = data\n",
" line.set_data(x,y)\n",
" return line,\n",
"\n",
"anim = animation.FuncAnimation(fig, animate, frames=un, interval=50)\n",
"display_animation(anim, default_mode='loop')"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"\n",
"\n",
"\n",
"
\n",
"\n",
"\n",
"\n"
],
"metadata": {},
"output_type": "pyout",
"prompt_number": 6,
"text": [
""
]
}
],
"prompt_number": 6
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##Lax Friedrichs Test 2\n",
"\n",
"Now we can run the same scheme but we'll change our CFL number so that $\\frac{\\Delta t}{\\Delta x}=0.5$"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"dt = .5*dx\n",
"u = u_ic() ##Reset our initial conditions (the square wave)\n",
"un = laxfriedrichs(u,nt,dt,dx)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 7
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we have our new data set, we can set up our next animation. We'll use all the same code from above, except now our `animate` function is already defined, so we can just set up a figure as before, then call `FuncAnimation`"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"fig = pyplot.figure();\n",
"ax = pyplot.axes(xlim=(0,4),ylim=(-.5,2));\n",
"line, = ax.plot([],[],lw=2);\n",
"\n",
"anim = animation.FuncAnimation(fig, animate, frames=un, interval=50)\n",
"display_animation(anim, default_mode='once')"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"\n",
"\n",
"\n",
"
\n",
"\n",
"\n",
"\n"
],
"metadata": {},
"output_type": "pyout",
"prompt_number": 18,
"text": [
""
]
}
],
"prompt_number": 18
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###Yikes.\n",
"\n",
"You may have noticed that we've made `nt` smaller than it has been for the other schemes. \n",
"With our CFL number equal to 0.5, this implicit method can handle about 60 iterations before it blows up to numbers that Python can't handle (which then get replaced with `Inf`, which then breaks the banded matrix solver).\n",
"\n",
"So we need to calm things down a little bit and prevent it from going nuclear so quickly.\n",
"\n",
"##Beam-Warming with damping"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Implicit methods can work well in many situations, but this isn't one of them. As you can see above, the undamped Beam-Warming quickly turns into nonsense. We can add damping to help keep things under control. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Add a damping term 'D' where:\n",
"\n",
"$$D = -\\epsilon_e (u^n_{i+2} - 4u^n_{i+1} + 6u^n_i - 4^n_{i-1} + u^n_{i-2})$$\n",
"\n",
"Note below that while we can easily include the damping term for most values of the `d` vector, `d[0]` and `d[1]` require a little extra attention. There are no periodic boundary conditions being used, so any value of $u_i$ where $i < 0$ is set equal to 1."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def dampit(u,eps,dt,dx):\n",
"\td = u[2]-.5*dt/dx*(u[3]**2/2-u[1]**2/2)+dt/(4*dx)*(u[3]**2-u[1]**2)\\\n",
"\t\t-eps*(u[4]-4*u[3]+6*u[2]-4*u[1]+u[0])\n",
"\treturn d\n",
"\n",
"def beamwarming_damp(u, nt, dt, dx):\n",
" ##Tridiagonal setup##\n",
" a = numpy.zeros_like(u)\n",
" b = numpy.ones_like(u)\n",
" c = numpy.zeros_like(u)\n",
" d = numpy.zeros_like(u)\n",
" \n",
" un = numpy.zeros((nt,len(u)))\n",
" un[:] = u.copy()\n",
" \n",
" eps = .125\n",
"\n",
" for n in range(1,nt): \n",
" u[0] = 1\n",
" E = utoE(u)\n",
" au = utoA(u)\n",
" \n",
" a[0] = -dt/(4*dx)*u[0]\n",
" a[1:] = -dt/(4*dx)*u[0:-1]\n",
" a[-1] = -dt/(4*dx)*u[-1]\n",
" \n",
" #b is all ones\n",
" \n",
" c[:-1] = dt/(4*dx)*u[1:]\n",
" \n",
" ###Calculate the damping factor for MOST of our u_vector\n",
" d[2:-2] = u[2:-2]-.5*dt/dx*(E[3:-1]-E[1:-3])+dt/(4*dx)\\\n",
" *(au[3:-1]-au[1:-3])\\\n",
" -eps*(u[4:]-4*u[3:-1]+6*u[2:-2]-4*u[1:-3]+u[:-4])\n",
" \n",
"\n",
" ###Calculate the damping factor for d[0] and d[1]\n",
" damp = numpy.concatenate((np.ones(2), u[:3]))\t\n",
" d[0] = dampit(damp,eps,dt,dx)\n",
" damp = numpy.concatenate((np.ones(1), u[:4]))\n",
" d[1] = dampit(damp,eps,dt,dx)\n",
" \n",
" ###subtract a[0]*LHS B.C to 'fix' thomas algorithm\n",
" d[0] = d[0] - u[0] * a[0] \n",
" \n",
" ab = numpy.matrix([c,b,a])\n",
" u = linalg.solve_banded((1,1), ab, d)\n",
" u[0]=1\n",
" un[n] = u.copy()\n",
" return un"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 19
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"u = u_ic()\n",
"sigma = 0.5\n",
"dt = sigma*dx\n",
"nt = 120\n",
"un = beamwarming_damp(u,nt,dt,dx)\n",
"\n",
"fig = pyplot.figure();\n",
"ax = pyplot.axes(xlim=(0,4),ylim=(-.5,2));\n",
"line, = ax.plot([],[],lw=2);\n",
"\n",
"anim = animation.FuncAnimation(fig, animate, frames=un, interval=50)\n",
"display_animation(anim, default_mode='once')"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"\n",
"\n",
"\n",
"