Script for solving advection problems in 1D using FDM

We will use two different schemes to solve the general advection problem in 1D

First we load some useful libraries

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from IPython.display import HTML
from matplotlib.patches import Rectangle

Next we define some physics constants

width   = 40           # width of domain
vel     = 4             # velocity in x direction
ampl    = 2             # amplitude of gaussian
sigma   = 1             # width of initial guassian
xc      = width/2.      # center of gaussian
dt      = 2e-2          # time step

and some numerical constants constants

nx      = 201                                          # Number of gridpoints in x-direction
nt      = 100                                          # Number of timesteps to compute
Xvec,dx = np.linspace(0, width, nx,  retstep=True)    # X coordinate vector, constant spacing

Set the initial conditions

C_init  = ampl*np.exp(-np.square(Xvec-xc)/sigma**2)
time    = 0

Make a function that computes temperature

# hide
# only store the latest time step solution to Cold
def advect1d_uw(Cold):
    Cnew=Cold.copy()
    if vel>0:
        for i in range(1,len(Xvec)):
            Cnew[i] =dt*(-vel* (Cold[i] - Cold[i-1])/dx) + Cold[i]
    else:
         for i in range(0,len(Xvec)-1):
            Cnew[i] =dt*(-vel* (Cold[i+1] - Cold[i])/dx) + Cold[i]

    return Cnew


def advect1d_ftcs(Cold):
    Cnew=Cold.copy()
    for i in range(1,len(Xvec)-1):
            Cnew[i] =dt*(-vel* (Cold[i+1] - Cold[i-1])/(2*dx)) + Cold[i]

    return Cnew

# write two functions that update the concentration using FTCS and upwind!

#def advect1d_uw(Cold):
#    Cnew=Cold.copy()
#    if vel>0:
#        for i in ???:
#            Cnew[i] ???
#    else:
#         for i in ???:
#            Cnew[i] =???
#
#    return Cnew


#def advect1d_ftcs(Cold):
#    Cnew=Cold.copy()
#    for i in ???:
#            Cnew[i] =???
#
#    return Cnew

Compute and visualize everything

# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure(figsize=(10,5))
ax = plt.axes(xlim=(0, width), ylim=(-0.2*ampl, 1.5*ampl))
ax.set_xlabel('Distance')
ax.set_ylabel('Concentration')
timeLabel=ax.text(0.02,0.98,'Time: ',transform=ax.transAxes,va='top')

num_solutions       = 3
lines = []
for i in range(0,num_solutions):
    line, = ax.plot([], [], lw=1)
    lines.append(line)

# Initialize Tnew
Cnew_uw=C_init
Cnew_ftcs=C_init

# Animation function which updates figure data.  This is called sequentially
def animate(i):
    timeLabel._text='Time: %.1f'%(i*dt)

    # use global keyword to store the latest solution and update it using fdm_solve function
    global Cnew_uw, Cnew_ftcs, Cnew_ana

    # numerical solutions
    if i==0:
        Cnew_uw=C_init
        Cnew_ftcs=C_init
        Cnew_ana=C_init
    else:
        Cnew_uw=advect1d_uw(Cnew_uw)
        Cnew_ftcs=advect1d_ftcs(Cnew_ftcs)
        Cnew_ana=ampl*np.exp(-np.square(Xvec-xc-i*dt*vel)/sigma**2)

    lines[0].set_data(Xvec, Cnew_uw)
    lines[1].set_data(Xvec, Cnew_ftcs)
    lines[2].set_data(Xvec, Cnew_ana)

    return lines


# Call the animator.
anim = animation.FuncAnimation(fig, animate,frames=nt,interval=50)
plt.close(anim._fig)

# Call function to display the animation
#HTML(anim.to_html5_video())  # lower resolution
HTML(anim.to_jshtml())  # higher resolution

Once Loop Reflect