{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Linear Advection with Finite Differences\n",
"## CH EN 6355 - Computational Fluid Dynamics\n",
"**Prof. Tony Saad (www.tsaad.net)
Department of Chemical Engineering
University of Utah**\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will be solving the linear advection equation\n",
"\\begin{equation}\n",
"\\frac{\\partial u}{\\partial t} = - c \\frac{\\partial u}{\\partial x}\n",
"\\end{equation}\n",
"subject to the initial condition $u(x,0)\\equiv u_0(x)$ on the domain $[0,L]$. The exact solution for this equation is given by\n",
"\\begin{equation}\n",
"u(x,t) = u_0(x-ct)\n",
"\\end{equation}\n"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib.animation as animation\n",
"plt.rcParams['animation.html'] = 'html5' # this is used to display animations in jupyter notebooks"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's start by defining an initial condition - let's take\n",
"\\begin{equation}\n",
" u_0(x) = \\sin(\\omega x)\n",
"\\end{equation}"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"n = 41\n",
"L = 1.0\n",
"dx = L/(n-1)\n",
"x = np.linspace(0,L,n)\n",
"u0 = np.sin(2.0*np.pi*x)\n",
"plt.plot(x,u0,'o')\n",
"plt.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Forward in Time, Central in Space (FTCS)\n",
"Using a forward Euler method for the time derivative and a central discretization for the spatial derivative, the FTCS scheme results in\n",
"\\begin{equation}\n",
"u_i^{n + 1} = u_i^n + c\\Delta t\\frac{u_{i + 1}^n - u_{i - 1}^n}{2\\Delta x}\n",
"\\end{equation}"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {},
"outputs": [],
"source": [
"# the following implementation does not use ghost cells\n",
"dt = 0.001 # s\n",
"tend = 2.0 # s\n",
"c = 10.0 # m/s - wave speed\n",
"cfl = c*dt/2.0/dx\n",
"\n",
"sol = []\n",
"sol.append(u0)\n",
"\n",
"t = 0.0\n",
"while t < tend:\n",
" un = sol[-1]\n",
" unew = un.copy()\n",
" unew[1:-1] = un[1:-1] - cfl * (un[2:] - un[:-2])\n",
" unew[-1] = un[-1] - cfl*(un[1] - un[-2]) # compute last point on the right using periodicity\n",
" unew[0] = unew[-1] # set periodic boundary on the left\n",
" sol.append(unew)\n",
" t += dt"
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 67,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"ims = []\n",
"fig = plt.figure(figsize=[4,3])\n",
"plt.grid()\n",
"\n",
"i = 0\n",
"for solution in sol:\n",
" if (i%10==0): # output frequency for frames \n",
" im = plt.plot(x,solution,'k-',animated=True)\n",
" plt.ylim(-1,1)\n",
" ims.append(im)\n",
" i+=1\n",
"ani = animation.ArtistAnimation(fig, ims, interval=35, blit=True,\n",
" repeat_delay=1000)\n",
"# ani.save('ftcs.mp4') \n",
"ani"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Forward in Time, Backward (Upwind) in Space (FTUS)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have shown in class that the previous FTCS scheme is unconditionally unstable. Here we will implement the upwind scheme in space"
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {},
"outputs": [],
"source": [
"tend = 1.0 # s\n",
"c = 1.0 # m/s - wave speed\n",
"dt = dx/c # run just at the CFL condition\n",
"cfl = c*dt/dx\n",
"\n",
"sol = []\n",
"sol.append(u0)\n",
"\n",
"t = 0.0\n",
"while t < tend:\n",
" un = sol[-1]\n",
" unew = un.copy()\n",
" unew[1:-1] = un[1:-1] - cfl * (un[1:-1] - un[:-2])\n",
" unew[-1] = un[-1] - cfl*(un[-1] - un[-2]) # compute last point on the right using periodicity\n",
" unew[0] = unew[-1] # set periodic boundary on the left\n",
" sol.append(unew)\n",
" t += dt"
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 72,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"ims = []\n",
"fig = plt.figure(figsize=[5,4],dpi = 150)\n",
"plt.grid()\n",
"\n",
"i = 0\n",
"for solution in sol:\n",
" if (i%2==0): \n",
" im = plt.plot(x,solution,'k-',animated=True)\n",
" plt.ylim(-1,1)\n",
" ims.append(im)\n",
" i+=1\n",
"ani = animation.ArtistAnimation(fig, ims, interval=35, blit=True,\n",
" repeat_delay=1000)\n",
"# ani.save('ftus.mp4') \n",
"ani"
]
},
{
"cell_type": "code",
"execution_count": 73,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"CSS style adapted from https://github.com/barbagroup/CFDPython. Copyright (c) Barba group\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 73,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import urllib\n",
"import requests\n",
"from IPython.core.display import HTML\n",
"def css_styling():\n",
" styles = requests.get(\"https://raw.githubusercontent.com/saadtony/NumericalMethods/master/styles/custom.css\")\n",
" return HTML(styles.text)\n",
"css_styling()\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python [default]",
"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.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}