Skip to content

Instantly share code, notes, and snippets.

@gforsyth

gforsyth/Burgers_Equation.ipynb

Last active October 27, 2015 21:39
Show Gist options
  • Save gforsyth/2062a86de03789120edd to your computer and use it in GitHub Desktop.
Save gforsyth/2062a86de03789120edd to your computer and use it in GitHub Desktop.
Download ZIP
Numerical Damping
Raw
README
Numerical Damping
Display the source blob
Display the rendered blob
Raw
Burgers_Equation.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Burgers Equation"
]
},
{
"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}\\left(\\frac{u^2}{2} \\right)$$\n",
"\n",
"We can also write this as:\n",
"\n",
"$$\\frac{\\partial u}{\\partial t} = - \\frac{\\partial F}{\\partial x}$$\n",
"\n",
"if we take $F = \\frac{u^2}{2}$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Initial Conditions\n",
"---\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"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy\n",
"from matplotlib import pyplot"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from matplotlib import animation\n",
"from JSAnimation.IPython_display import display_animation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##### Task 1:\n",
"\n",
"Complete the function `u_initial` so that it returns a NumPy array with the initial conditions specified above."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def u_initial():\n",
" \n",
" \n",
" return u"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"nx = 81\n",
"nt = 70\n",
"dx = 4.0/(nx-1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"computeF = lambda u: (u/2)**2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##### Task 2:\n",
"\n",
"What does `computeF` do? "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## MacCormack Scheme\n",
"\n",
"$$u^*_i = u^n_i - \\frac{\\Delta t}{\\Delta x} (F^n_{i+1}-F^n_{i}) \\ \\ \\ \\ \\ \\ (predictor)$$\n",
"\n",
"$$u^{n+1}_i = \\frac{1}{2} \\left(u^n_i + u^*_i - \\frac{\\Delta t}{\\Delta x} (F^*_i - F^{*}_{i-1}) \\right) \\ \\ \\ \\ \\ \\ (corrector)$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##### Task 3: \n",
"\n",
"Complete the `maccormack` function below using array operations"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def maccormack(u, nt, dt, dx):\n",
" un = numpy.zeros((nt,len(u)))\n",
" un[:] = u.copy()\n",
" ustar = u.copy()\n",
" \n",
" for i in range(1,nt):\n",
" F = computeF(u)\n",
" \n",
" ustar = \n",
" \n",
" Fstar = computeF(ustar)\n",
" \n",
" un[i] = \n",
" \n",
" u = un[i].copy()\n",
" \n",
" return un"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### CFL = 1"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def animate(data):\n",
" x = numpy.linspace(0,4,nx)\n",
" y = data\n",
" line.set_data(x,y)\n",
" return line,"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"u = u_initial()\n",
"sigma = 1\n",
"dt = sigma*dx\n",
"\n",
"un = maccormack(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')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Numerical damping"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The oscillations at the front of the shock are characteristic of all 2nd-order methods, but there are ways to remove them. One option is numerical damping. \n",
"\n",
"A common damping term used with MacCormack is\n",
"\n",
"$$\\epsilon \\left(u^n_{i+1} - 2u^n_i + u^n_{i-1} \\right)$$\n",
"\n",
"##### Task 4\n",
"\n",
"Add this term to the *predictor* step of your MacCormack function above and experiment with values of $\\epsilon$ from $0 < \\epsilon < 1$. Try to find a value of $\\epsilon$ that damps out the oscillations without messing up the rest of the solution too much. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
"<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\n",
"<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\n",
"<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>\n",
"<link href='http://fonts.googleapis.com/css?family=Nixie+One' rel='stylesheet' type='text/css'>\n",
"<style>\n",
"\n",
"@font-face {\n",
" font-family: \"Computer Modern\";\n",
" src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
"}\n",
"\n",
"#notebook_panel { /* main background */\n",
" background: rgb(245,245,245);\n",
"}\n",
"\n",
"div.cell { /* set cell width */\n",
" width: 750px;\n",
"}\n",
"\n",
"div #notebook { /* centre the content */\n",
" background: #fff; /* white background for content */\n",
" width: 1000px;\n",
" margin: auto;\n",
" padding-left: 0em;\n",
"}\n",
"\n",
"#notebook li { /* More space between bullet points */\n",
"margin-top:0.8em;\n",
"}\n",
"\n",
"/* draw border around running cells */\n",
"div.cell.border-box-sizing.code_cell.running { \n",
" border: 1px solid #111;\n",
"}\n",
"\n",
"/* Put a solid color box around each cell and its output, visually linking them*/\n",
"div.cell.code_cell {\n",
" background-color: rgb(256,256,256); \n",
" border-radius: 0px; \n",
" padding: 0.5em;\n",
" margin-left:1em;\n",
" margin-top: 1em;\n",
"}\n",
"\n",
"div.text_cell_render{\n",
" font-family: 'Alegreya Sans' sans-serif;\n",
" line-height: 140%;\n",
" font-size: 125%;\n",
" font-weight: 400;\n",
" width:600px;\n",
" margin-left:auto;\n",
" margin-right:auto;\n",
"}\n",
"\n",
"\n",
"/* Formatting for header cells */\n",
".text_cell_render h1 {\n",
" font-family: 'Nixie One', serif;\n",
" font-style:regular;\n",
" font-weight: 400; \n",
" font-size: 45pt;\n",
" line-height: 100%;\n",
" color: rgb(0,51,102);\n",
" margin-bottom: 0.5em;\n",
" margin-top: 0.5em;\n",
" display: block;\n",
"}\t\n",
".text_cell_render h2 {\n",
" font-family: 'Nixie One', serif;\n",
" font-weight: 400;\n",
" font-size: 30pt;\n",
" line-height: 100%;\n",
" color: rgb(0,51,102);\n",
" margin-bottom: 0.1em;\n",
" margin-top: 0.3em;\n",
" display: block;\n",
"}\t\n",
"\n",
".text_cell_render h3 {\n",
" font-family: 'Nixie One', serif;\n",
" margin-top:16px;\n",
"\tfont-size: 22pt;\n",
" font-weight: 600;\n",
" margin-bottom: 3px;\n",
" font-style: regular;\n",
" color: rgb(102,102,0);\n",
"}\n",
"\n",
".text_cell_render h4 { /*Use this for captions*/\n",
" font-family: 'Nixie One', serif;\n",
" font-size: 14pt;\n",
" text-align: center;\n",
" margin-top: 0em;\n",
" margin-bottom: 2em;\n",
" font-style: regular;\n",
"}\n",
"\n",
".text_cell_render h5 { /*Use this for small titles*/\n",
" font-family: 'Nixie One', sans-serif;\n",
" font-weight: 400;\n",
" font-size: 16pt;\n",
" color: rgb(163,0,0);\n",
" font-style: italic;\n",
" margin-bottom: .1em;\n",
" margin-top: 0.8em;\n",
" display: block;\n",
"}\n",
"\n",
".text_cell_render h6 { /*use this for copyright note*/\n",
" font-family: 'PT Mono', sans-serif;\n",
" font-weight: 300;\n",
" font-size: 9pt;\n",
" line-height: 100%;\n",
" color: grey;\n",
" margin-bottom: 1px;\n",
" margin-top: 1px;\n",
"}\n",
"\n",
".CodeMirror{\n",
" font-family: \"PT Mono\";\n",
" font-size: 90%;\n",
"}\n",
"\n",
"</style>\n",
"<script>\n",
" MathJax.Hub.Config({\n",
" TeX: {\n",
" extensions: [\"AMSmath.js\"],\n",
" equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
" },\n",
" tex2jax: {\n",
" inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
" displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
" },\n",
" displayAlign: 'center', // Change this to 'center' to center equations.\n",
" \"HTML-CSS\": {\n",
" styles: {'.MathJax_Display': {\"margin\": 4}}\n",
" }\n",
" });\n",
"</script>\n"
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.core.display import HTML\n",
"def css_styling():\n",
" styles = open(\"numericalmoocstyle.css\", \"r\").read()\n",
" return HTML(styles)\n",
"css_styling()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.4.3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Raw
numericalmoocstyle.css
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Nixie+One' rel='stylesheet' type='text/css'>
<style>
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}
#notebook_panel { /* main background */
background: rgb(245,245,245);
}
div.cell { /* set cell width */
width: 750px;
}
div #notebook { /* centre the content */
background: #fff; /* white background for content */
width: 1000px;
margin: auto;
padding-left: 0em;
}
#notebook li { /* More space between bullet points */
margin-top:0.8em;
}
/* draw border around running cells */
div.cell.border-box-sizing.code_cell.running {
border: 1px solid #111;
}
/* Put a solid color box around each cell and its output, visually linking them*/
div.cell.code_cell {
background-color: rgb(256,256,256);
border-radius: 0px;
padding: 0.5em;
margin-left:1em;
margin-top: 1em;
}
div.text_cell_render{
font-family: 'Alegreya Sans' sans-serif;
line-height: 140%;
font-size: 125%;
font-weight: 400;
width:600px;
margin-left:auto;
margin-right:auto;
}
/* Formatting for header cells */
.text_cell_render h1 {
font-family: 'Nixie One', serif;
font-style:regular;
font-weight: 400;
font-size: 45pt;
line-height: 100%;
color: rgb(0,51,102);
margin-bottom: 0.5em;
margin-top: 0.5em;
display: block;
}
.text_cell_render h2 {
font-family: 'Nixie One', serif;
font-weight: 400;
font-size: 30pt;
line-height: 100%;
color: rgb(0,51,102);
margin-bottom: 0.1em;
margin-top: 0.3em;
display: block;
}
.text_cell_render h3 {
font-family: 'Nixie One', serif;
margin-top:16px;
font-size: 22pt;
font-weight: 600;
margin-bottom: 3px;
font-style: regular;
color: rgb(102,102,0);
}
.text_cell_render h4 { /*Use this for captions*/
font-family: 'Nixie One', serif;
font-size: 14pt;
text-align: center;
margin-top: 0em;
margin-bottom: 2em;
font-style: regular;
}
.text_cell_render h5 { /*Use this for small titles*/
font-family: 'Nixie One', sans-serif;
font-weight: 400;
font-size: 16pt;
color: rgb(163,0,0);
font-style: italic;
margin-bottom: .1em;
margin-top: 0.8em;
display: block;
}
.text_cell_render h6 { /*use this for copyright note*/
font-family: 'PT Mono', sans-serif;
font-weight: 300;
font-size: 9pt;
line-height: 100%;
color: grey;
margin-bottom: 1px;
margin-top: 1px;
}
.CodeMirror{
font-family: "PT Mono";
font-size: 90%;
}
</style>
<script>
MathJax.Hub.Config({
TeX: {
extensions: ["AMSmath.js"],
equationNumbers: { autoNumber: "AMS", useLabelIds: true}
},
tex2jax: {
inlineMath: [ ['$','$'], ["\\(","\\)"] ],
displayMath: [ ['$$','$$'], ["\\[","\\]"] ]
},
displayAlign: 'center', // Change this to 'center' to center equations.
"HTML-CSS": {
styles: {'.MathJax_Display': {"margin": 4}}
}
});
</script>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment