{
"cells": [
{
"cell_type": "raw",
"metadata": {},
"source": [
"Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved BSD-3 license. (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[@LorenaABarba](https://twitter.com/LorenaABarba)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"12 steps to Navier–Stokes\n",
"=====\n",
"***"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This lesson complements the first interactive module of the online [CFD Python](https://github.com/barbagroup/CFDPython) class, by Prof. Lorena A. Barba, called **12 Steps to Navier–Stokes.** It was written with BU graduate student Gilbert Forsyth."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Array Operations with NumPy\n",
"----------------\n",
"\n",
"For more computationally intensive programs, the use of built-in Numpy functions can provide an increase in execution speed many-times over. As a simple example, consider the following equation:\n",
"\n",
"$$u^{n+1}_i = u^n_i-u^n_{i-1}$$\n",
"\n",
"Now, given a vector $u^n = [0, 1, 2, 3, 4, 5]\\ \\ $ we can calculate the values of $u^{n+1}$ by iterating over the values of $u^n$ with a for loop. "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1\n",
"1\n",
"1\n",
"1\n",
"1\n"
]
}
],
"source": [
"u = numpy.array((0, 1, 2, 3, 4, 5))\n",
"\n",
"for i in range(1, len(u)):\n",
" print(u[i] - u[i-1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is the expected result and the execution time was nearly instantaneous. If we perform the same operation as an array operation, then rather than calculate $u^n_i-u^n_{i-1}\\ $ 5 separate times, we can slice the $u$ array and calculate each operation with one command:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1, 1, 1, 1, 1])"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u[1:] - u[0:-1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What this command says is subtract the 0th, 1st, 2nd, 3rd, 4th and 5th elements of $u$ from the 1st, 2nd, 3rd, 4th, 5th and 6th elements of $u$. \n",
"\n",
"### Speed Increases\n",
"\n",
"For a 6 element array, the benefits of array operations are pretty slim. There will be no appreciable difference in execution time because there are so few operations taking place. But if we revisit 2D linear convection, we can see some substantial speed increases. \n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"nx = 81\n",
"ny = 81\n",
"nt = 100\n",
"c = 1\n",
"dx = 2 / (nx - 1)\n",
"dy = 2 / (ny - 1)\n",
"sigma = .2\n",
"dt = sigma * dx\n",
"\n",
"x = numpy.linspace(0, 2, nx)\n",
"y = numpy.linspace(0, 2, ny)\n",
"\n",
"u = numpy.ones((ny, nx)) ##create a 1xn vector of 1's\n",
"un = numpy.ones((ny, nx)) \n",
"\n",
"###Assign initial conditions\n",
"\n",
"u[int(.5 / dy): int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With our initial conditions all set up, let's first try running our original nested loop code, making use of the iPython \"magic\" function `%%timeit`, which will help us evaluate the performance of our code. \n",
"\n",
"**Note**: The `%%timeit` magic function will run the code several times and then give an average execution time as a result. If you have any figures being plotted within a cell where you run `%%timeit`, it will plot those figures repeatedly which can be a bit messy. \n",
"\n",
"The execution times below will vary from machine to machine. Don't expect your times to match these times, but you _should_ expect to see the same general trend in decreasing execution time as we switch to array operations."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3.07 s ± 15.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
]
}
],
"source": [
"%%timeit\n",
"u = numpy.ones((ny, nx))\n",
"u[int(.5 / dy): int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2\n",
"\n",
"for n in range(nt + 1): ##loop across number of time steps\n",
" un = u.copy()\n",
" row, col = u.shape\n",
" for j in range(1, row):\n",
" for i in range(1, col):\n",
" u[j, i] = (un[j, i] - (c * dt / dx * \n",
" (un[j, i] - un[j, i - 1])) - \n",
" (c * dt / dy * \n",
" (un[j, i] - un[j - 1, i])))\n",
" u[0, :] = 1\n",
" u[-1, :] = 1\n",
" u[:, 0] = 1\n",
" u[:, -1] = 1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the \"raw\" Python code above, the mean execution time achieved was 3.07 seconds (on a MacBook Pro Mid 2012). Keep in mind that with these three nested loops, that the statements inside the **j** loop are being evaluated more than 650,000 times. Let's compare that with the performance of the same code implemented with array operations:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"7.38 ms ± 105 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
]
}
],
"source": [
"%%timeit\n",
"u = numpy.ones((ny, nx))\n",
"u[int(.5 / dy): int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2\n",
"\n",
"for n in range(nt + 1): ##loop across number of time steps\n",
" un = u.copy()\n",
" u[1:, 1:] = (un[1:, 1:] - (c * dt / dx * (un[1:, 1:] - un[1:, 0:-1])) -\n",
" (c * dt / dy * (un[1:, 1:] - un[0:-1, 1:])))\n",
" u[0, :] = 1\n",
" u[-1, :] = 1\n",
" u[:, 0] = 1\n",
" u[:, -1] = 1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As you can see, the speed increase is substantial. The same calculation goes from 3.07 seconds to 7.38 milliseconds. 3 seconds isn't a huge amount of time to wait, but these speed gains will increase exponentially with the size and complexity of the problem being evaluated. "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.core.display import HTML\n",
"def css_styling():\n",
" styles = open(\"../styles/custom.css\", \"r\").read()\n",
" return HTML(styles)\n",
"css_styling()"
]
}
],
"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.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 1
}