For this comparison I'm going to assume most of you are primarily Matlab users. Matlab is great, especially in a university environment. It's an easy to use, interpreted, high-level language with automatic memory management and lots of supporting libraries. Python shares those same advantages. What's the difference then? Python has some benefits that may become more important to you as your problems get larger, or as you move into industry. Some of these benefits include:
Are there any drawbacks? The only one I've come across is Simulink. I don't know of anything of comparable capability in the Python world. But, since this is a fluid dynamics conference, probably none of us is worried about that.
Some packages you may want to install include (you will definitely need the first three):
Many exist, and I'm not familiar enough with them all to recommend one over another. I like the one in the official Python docs. It's probably more detailed than you want for a first exposure, but it is a good resource.
Here are two useful resources called NumPy for Matlab Users: one, two.
Today we will go through two simple examples to introduce you to the syntax and to show how to wrap an external Fortran/C code. If you are viewing this in nbviewer I suggest you download it first (button in upper right corner), then open it up at https://try.jupyter.org so you can edit and follow along. Later you can install IPython so you can make edits locally rather than on a remote server. First, evaluate the cell below. Press shift-enter to evaluate a cell. This cell imports some libraries we will need. The first line is only for IPython (you wouldn't use it in a normal Python script). It just tells IPython that we want to see plots inline. The next line imports numpy, which contains a lot of important array operations. The last imports matplotlib, which is a plotting package.
In [2]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
We are going to compute induced and parasite drag (in a very basic way) as an example to introduce scripting. \begin{align} q&= \frac{1}{2} \rho V_\infty^2 \\ D_i &= \frac{L^2}{q \pi b^2 e} \\ D_p &= {C_D}_p q S \end{align}
The first few lines make some imports. The first line is only for Jupyter notebooks. It just tells the notebook to show plots inline. You wouldn't need that in a normal python script.
In [3]:
from math import pi
# all in standard English units
# atmosphere
rho = 0.0024 # air density
# geometry
b = 8.0 # wing span
chord = 1.0
# mass properties
W = 2.4 # total weight of aircraft
# other parameters
e = 0.9 # Oswald efficiency factor
CDp = 0.02 # could compute but just input for simplicity
# an array of wind speeds
V = np.linspace(10, 30, 100)
# Induced drag
q = 0.5*rho*V**2 # dyamic pressure
L = W # equilibrium flight
Di = L**2/(q*pi*b**2*e)
# parasite drag
S = b*chord
Dp = CDp*q*S
# these next 3 lines purely for style in the plot (loading a predefined styleshet)
# I have my own custom styles I use, but for this example let's use one of matplotlib's
plt.style.use('ggplot')
plt.rcParams.update({'font.size': 16})
colors = plt.rcParams['axes.color_cycle'] # grab the current color scheme
# plot it
plt.figure()
plt.plot(V, Di)
plt.plot(V, Dp)
plt.plot(V, Di+Dp)
plt.xlabel('V (ft/s)')
plt.ylabel('Drag (lbs)')
# label the plots
plt.text(25, 0.06, 'induced drag', color=colors[0])
plt.text(12, 0.06, 'parasite drag', color=colors[1])
plt.text(20, 0.17, 'total drag', color=colors[2])
Out[3]:
Try it yourself. Let's do the same calculation, but with reusable functions. In Python functions are easy to define. A simple example is below (note that unlike Matlab, you can have as many functions as you want in a file.)
In [4]:
def func(x, y):
add = x + y
mult = x * y
return add, mult
a, m = func(1.0, 3.0)
print 'a =', a, 'm =', m
a, m = func(2.0, 7.0)
print 'a =', a, 'm =', m
In [5]:
def induced_drag():
pass
def parasite_drag():
pass
# atmosphere
rho = 0.0024 # air density
# geometry
b = 8.0 # wing span
chord = 1.0
# mass properties
W = 2.4 # total weight of aircraft
# other parameters
e = 0.9 # Oswald efficiency factor
CDp = 0.02 # could compute but just input for simplicity
# wind speeds
V = np.linspace(10, 30, 100)
Finally, let's do it once more, but in an object-oriented style.
In [6]:
class UAV(object):
def __init__(self, b, chord, W, rho):
self.b = b
self.S = b*chord
self.L = W
self.rho = rho
def induced_drag(self, V, e):
q = 0.5*self.rho*V**2
Di = self.L**2/(q*pi*self.b**2*e)
return Di
def parasite_drag(self, V, CDp):
q = 0.5*self.rho*V**2
Dp = CDp*q*self.S
return Dp
# atmosphere
rho = 0.0024 # air density
# geometry
b = 8.0 # wing span
chord = 1.0
# mass properties
W = 2.4 # total weight of aircraft
# setup UAV object
uav = UAV(b, chord, W, rho)
# setup sweep
V = np.linspace(10, 30, 100)
# idrag
e = 0.9 # Oswald efficiency factor
Di = uav.induced_drag(V, e)
# pdrag
CDp = 0.02
Dp = uav.parasite_drag(V, CDp)
# style
plt.style.use('fivethirtyeight')
# plot it
plt.figure()
plt.plot(V, Di)
plt.plot(V, Dp)
plt.plot(V, Di+Dp)
plt.xlabel('V (ft/s)')
plt.ylabel('Drag (lbs)')
# label the plots
colors = plt.rcParams['axes.color_cycle']
plt.text(25, 0.06, 'induced drag', color=colors[0])
plt.text(12, 0.06, 'parasite drag', color=colors[1])
plt.text(20, 0.17, 'total drag', color=colors[2])
Out[6]:
This next example is going to solve Laplace's equation on a grid. First, we will do it in pure Python, then we will rewrite a portion of the code in Fortran and call it in Python for improved speed. Recall Laplace's equation:
$$ \nabla^2 \phi = 0 $$where $\phi is some scalar. For a regular rectangular grid, with equal spacing in x and y, you might recall that a simple iterative method for solving this equation consists of the following update rule:
$$ \phi_{i, j} = \frac{1}{4} (\phi_{i+1, j} + \phi_{i-1, j} + \phi_{i, j+1} + \phi_{i, j-1})$$In other words, each cell updates its value using the average value of all of its neighbors (note that they are much more efficient ways to solve Laplace's equation on a grid. For our purpose we just want to keep things simple). This process must be repeated for every cell in the domain, and repeated until converged.
We are going to run a simple case where boundary values are provided at the top, bottom, left, and right edges. You should iterate until the maximum change in $\phi$ is below some tolerance (tol) or until a maximum number of iterations is reached (iter_max). n is the number of cells (same discretization in x and y).
I've started a script for you below. See if you can fill in the details. I've not provided all the syntax you will need to know so you may have to look some things up. A full implementation is down below, but don't peek unless you are really stuck!
In [ ]:
from math import fabs
def laplace_grid_python(n, top, bottom, left, right, tol, iter_max):
# initialize
phi = np.zeros((n+1, n+1))
iters = 0 # number of iterations
err_max = 1e6 # maximum error in grid (start at some arbitrary number just to enter loop)
# set boundary conditions
# run while loop until tolerance reached or max iterations
while ():
# reset the maximum error to something small (I suggest something like -1)
err_max = -1.0
# loop over all *interior* cells
for i in range():
for j in range():
# save previous point for computing error later
phi_prev = phi[i, j]
# update point
phi[i, j] =
# update maximum error
err_max =
# update iteration count
iters += 1
return phi, err_max, iters
# run a sample case (50 x 50 grid with bottom and left and 1.0, top and right at 0.0)
n = 50
top = 0.0
bottom = 1.0
left = 1.0
right = 0.0
tol = 1e-5
iter_max = 10000
phi, err_max, iters = laplace_grid_python(n, top, bottom, left, right, tol, iter_max)
# plot it
x = np.linspace(0, 1, n+1)
y = np.linspace(0, 1, n+1)
[X, Y] = np.meshgrid(x, y)
plt.figure()
plt.contourf(X, Y, phi, 100, cmap=plt.cm.get_cmap("YlGnBu"))
plt.colorbar()
plt.show()
In our implementation below we are adding the IPython flag %%timeit at the top. This will run the whole block some number of times and report the best time back.
Adding some blank space below just to visually separate the answer.
In [10]:
%%timeit
from math import fabs
def laplace_grid_python(n, top, bottom, left, right, tol, iter_max):
# initialize
phi = np.zeros((n+1, n+1))
iters = 0 # number of iterations
err_max = 1e6 # maximum error in grid (start at some arbitrary number just to enter loop)
# set boundary conditions
phi[0, :] = bottom
phi[-1, :] = top
phi[:, 0] = left
phi[:, -1] = right
# run while loop until tolerance reached or max iterations
while (err_max > tol and iters < iter_max):
# reset the maximum error to something small (I suggest something like -1)
err_max = -1.0
# loop over all interior cells
for i in range(1, n):
for j in range(1, n):
# save previous point
phi_prev = phi[i, j]
# update point
phi[i, j] = (phi[i-1,j] + phi[i+1,j] + phi[i,j-1] + phi[i,j+1])/4.0
# update maximum error
err_max = max(err_max, fabs(phi[i, j] - phi_prev))
# update iteration count
iters += 1
return phi, err_max, iters
# run a sample case (50 x 50 grid with bottom and left and 1.0, top and right at 0.0)
n = 50
top = 0.0
bottom = 1.0
left = 1.0
right = 0.0
tol = 1e-5
iter_max = 10000
phi, err_max, iters = laplace_grid_python(n, top, bottom, left, right, tol, iter_max)
# plot it
x = np.linspace(0, 1, n+1)
y = np.linspace(0, 1, n+1)
[X, Y] = np.meshgrid(x, y)
plt.figure()
plt.contourf(X, Y, phi, 100, cmap=plt.cm.get_cmap("YlGnBu"))
plt.colorbar()
plt.show()
This takes a while. Let's move the double for loop computation to Fortran. I've supplied a file called laplace.f90 where I've done this for you. We just need to build this as a shared library so we can call it from Python. Open a terminal (you can do this in try.jupyter.org as well). We will compile the fortran code to a shared library with f2py, (can also use standard gfortran compilation but you get more flexibility and automatic setup with f2py). In all of the below I am using an O2 optimization flag. Note the underscore in the shared library name. This is just convention. Note that if you make a mistake in importing, IPython caches your modules so you'd need to restart the Kernel. You can even do all of this <try.jupyter.org> by opening a terminal.
Using f2py
f2py -c --opt=-O2 -m _laplace laplace.f90
Using setup script. Open up a file and call it setup.py. At a minimum this is all it needs:
from numpy.distutils.core import setup, Extension
setup(
ext_modules=[Extension('_laplace', ['laplace.f90'], extra_compile_args=['-O2'])]
)
Usually a lot more information would be added (name, license, other python packages, etc.) You can read more about setuptools later.
To build it one normally just uses build or install commands, but we will built it inplace for testing.
python setup.py build_ext --inplace
Now we can call it from Python just as a regular method. An example is shown below doing the exact same thing as before, but calling the Fortran code in laplacegridfortran. This runs over 10x faster (and could be even faster if we used ifort instead of gfortran). How much faster the code is of course depends ont he problem as you change n the difference will either become more or less significant.
In [11]:
%%timeit
from _laplace import laplacegridfortran
n = 50
top = 0.0
bottom = 1.0
left = 1.0
right = 0.0
tol = 1e-5
iter_max = 10000
phi, err_max, iters = laplacegridfortran(n, top, bottom, left, right, tol, iter_max)
# plot it
x = np.linspace(0, 1, n+1)
y = np.linspace(0, 1, n+1)
[X, Y] = np.meshgrid(x, y)
plt.figure()
plt.contourf(X, Y, phi, 100, cmap=plt.cm.get_cmap("YlGnBu"))
plt.colorbar()
plt.show()
In [ ]: