Step 1: Import packages

  • Numpy` is used for the vectors
  • scipy is unused here, but is common
  • matpllotlib.pyplot is used for plotting
In [51]:
%matplotlib inline
import numpy as np
import scipy as sp
import scipy.sparse as sparse
import matplotlib.pyplot as plt
import matplotlib.cm as cm

Step 2: Problem Setup

In [88]:
n = 200
c = 1.0
L = 1.0
dx = L / n
dxi = 1.0 / dx

Ep = sparse.diags([1,1],[1,1-n],shape=(n, n), format='csr')
Em = Ep.T
Du = dxi * (Ep - Em) / 2 # Centered, 2nd-order upwind, or 3rd-order upwind?

Tfinal = 1.0
CFL = 0.13
dt = CFL * (dx / c)
Nsteps = np.ceil(Tfinal / dt)
dt = Tfinal / Nsteps

Step 3: Initial Conditions

In [79]:
def u0_sin(x):
    return np.sin(4*np.pi*x)

def u0_gaussian_pulse(x):
    return np.sin(np.pi*x)**50

def u0_exp(x):
    argx = -1.0*(x - 0.5)**2 / 0.001
    return np.exp(argx)

Step 4: Choose Initial Condition

In [108]:
x = dx*np.linspace(1,n,n)

# Uncomment an initial condition to use it
#u = u0_sin(x)
#u = u0_gaussian_pulse(x)
u = u0_exp(x)

Step 5: Setup Solver and Plot

In [109]:
plt.figure(figsize=(10,8))
plt.plot(x, u, label='initial', color='black', marker='+', linestyle=':',\
         linewidth=2.5)

f1 = 0*u
f2 = 0*u

iostep = np.floor(Nsteps/10.0)

colors = cm.rainbow(np.linspace(0, 1, 11))
color_ctr = 0

for k in range(int(Nsteps)):
    if k==0:
        c0 = 1.0
        c1 = 0.0
        c2 = 0.0
    elif k==1:
        c0 = 1.5
        c1 = -0.5
        c2 = 0.0
    else:
        c0 = 23.0/12.0
        c1 = -16.0/12.0
        c2 = 5.0/12.0
    
    f0 = -Du*u
    f = c0*f0 + c1*f1 + c2*f2
    f2 = f1
    f1 = f0
    u = u + dt*f
    
    if k%iostep == 0:
        plt.plot(x, u, label='timestep '+str(k), color=colors[color_ctr],\
                 linestyle='--', linewidth=2.5)
        color_ctr += 1
        
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))

#Uncomment the appropriate y limits for your initial condition
#plt.ylim(-1.2, 1.2) # sin initial condition ylim
#plt.ylim(-0.2, 1.2) # gaussian pulse ylim
plt.ylim(-0.5, 1.2) # exp ylim

plt.xlim(0.0, 1.0)
plt.xlabel("--x--")
plt.ylabel("u(x,t)")
plt.title('1D Advection: 2nd-Order Finite difference, N=200')
Out[109]:
<matplotlib.text.Text at 0x7f402b85fd90>
In [ ]: