In [1]:
%matplotlib inline
In [2]:
from pylab import *
In [3]:
u = 2.
In [4]:
xlower = 0.
xupper = 1.
num_cells = 20
dx = (xupper-xlower)/num_cells
x = arange(xlower - dx/2, xupper + dx, dx)
In [5]:
dt = 0.5 * dx/u
cfl = u*dt/dx
print(cfl)
0.5
In [6]:
def upwind_step(Qn):
Qnp = Qn.copy()
# periodic BCs:
Qn[0] = Qn[-2]
Qn[-1] = Qn[1]
for i in range(1, len(x)-1):
Qnp[i] = Qn[i] - cfl*(Qn[i] - Qn[i-1])
return Qnp
In [7]:
Qn = sin(2*pi*x)
plot(x[1:-1], Qn[1:-1], 'bo-')
grid(True)
xlim(xlower, xupper)
Out[7]:
(0.0, 1.0)
In [8]:
def plotQ(Qn, tn):
plot(x[1:-1], Qn[1:-1], 'bo-')
grid(True)
xlim(xlower, xupper)
title('Time t = %.3f' % tn)
In [9]:
tn = 0.
Qn = where(x<0.5, 1., 0.)
plotQ(Qn,tn)
In [10]:
tn = tn+dt
Qn = upwind_step(Qn)
figure(figsize=(6,3))
plotQ(Qn,tn)
In [11]:
#Qn = sin(2*pi*x)
Qn = where(x<0.5, 1., 0.)
tn = 0.
figure(figsize=(6,3))
plotQ(Qn,tn)
for n in range(5):
tn = tn + dt
Qn = upwind_step(Qn)
plotQ(Qn,tn)
Note on numpy arrays¶
The bug in my code yesterday was that I set
Qn = where(x<0.5, 1, 0)
rather than
Qn = where(x<0.5, 1., 0.)
The problem is that using integers 1, 0 rather than floats 1, 0. results in Qn being a numpy array with dtype=int (an array of integers). Then if you modify one element, e.g.
Qn[0] = 0.6
since the array is hardwared to only hold integers, this gets rounded up to 1 instead.
In [ ]: