The classical advection equation is given by: $$ \begin{equation} \label{eq:6101} \frac{\partial u}{\partial t} + a_0 \frac{\partial u}{\partial x} = 0, \end{equation} $$
where \( u = u(x, t) \). It is linear and also have an analytical solution: $$ \begin{equation}\label{eq:hyp_analytical} u = u_0(x-a_0\,t ), \end{equation} $$
In which \( u_0(x) \) is an initial solution/wave (\( u(x, 0)=u_0(x) \)), which we choose to be a smooth sine squared wave between \( x=0 \) and \( x=0.2 \): $$ \begin{align}\label{eq:initialWave} u_0(x) = sin^2(\pi(x/0.2)), \quad 0 < x < 0.2, \\ u_0(x) = 0, \quad 0.2 \leq x \leq l, \label{_auto1} \end{align} $$
where \( l = 1 \). Choose a suitable discretization \( \Delta x = \frac{l}{N} \) and calculate \( \Delta t \) from the CFL condition. Simulate from \( t_{start}=0 \) to \( t_{end}=0.8 \). This way we avoid the need to introduce special boundary conditions. (The wave will not leave the domain, and \( u(0)=0 \) and \( u(l)=0 \)). Otherwise the following boundarycondtions should be applied at the right boundary: $$ \begin{equation}\label{eq:bc} u(l)^{n + 1} = u(l - a \cdot \Delta t)^n, \end{equation} $$
where \( n \) denots the previous timestep, \( a \) is the wavespeed (\( a_0 \)), and \( u(l - a \cdot \Delta t) \) can be found by interpolation.
The forward in time backward in space (ftbs), or upwind scheme in non-conservative form may be written as follows: $$ \begin{equation} \label{eq:ftbs_conservative} u_j^{n+1}=u_j^n - \frac{a_0 \Delta t}{\Delta x}\left(u_j^n-u_{j-1}^n\right) \end{equation} $$
a) Compare ftbs with the analytical solution for the linear advection equation. choose \( a_0=1 \) and a cfl-constraint condition of \( c=1 \).
Now let's look at a nonlinear case of the advection equation, in which the wavespeed a is nolonger constant: $$ \begin{equation} \label{eq:advection_nonLin} \frac{\partial u}{\partial t} + a(u) \frac{\partial u}{\partial x} = \frac{\partial u}{\partial t} + (0.9 + 0.1u) \frac{\partial u}{\partial x} = 0 \end{equation} $$
Now we can modify the nonconservative ftbs scheme to work for the non-constant wavespeed: $$ \begin{equation} \label{eq:ftbs_conservative_nonLin} u_j^{n+1}=u_j^n - \frac{a_j^n \Delta t}{\Delta x}\left(u_j^n-u_{j-1}^n\right) \end{equation} $$
b) set \( a(u) = 0.9 + 0.1u \) and solve Eq. \eqref{eq:advection_nonLin} with the scheme given in \eqref{eq:ftbs_conservative_nonLin}
c) Rewrite Eq. \eqref{eq:advection_nonLin} so that it is written in conservative form: $$ \begin{equation} \label{eq:advection_nonLin2} \frac{\partial u}{\partial t} + \frac{\partial F}{\partial x} = 0 \end{equation} $$
What is the expression of \( F \) if \( a(u) = 0.9 + 0.1u \)? Solve Eq. \eqref{eq:advection_nonLin2} numerically with a conservative ftbs scheme: $$ \begin{equation} \label{eq:ftbs_conservativeFlux} u_j^{n+1}=u_j^n - \frac{\Delta t}{\Delta x}\left(F_j^n-F_{j-1}^n\right) \end{equation} $$
d) Do a), b), and c) but using the MacCormack scheme instead of ftbs
# python_exercises/pexercise6/problem1.py
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from scipy import interpolate
LNWDT=2; FNT=15
plt.rcParams['lines.linewidth'] = LNWDT; plt.rcParams['font.size'] = FNT
# function defining the initial condition
# def f(x):
# """Assigning a value of 1.0 for values less than 0.1"""
# f = np.zeros_like(x)
# f[np.where(x <= 0.1)] = 1.0
# return f
def f(x):
"""A smooth sin^2 function between x_left and x_right"""
f = np.zeros_like(x)
x_left = 0.0
x_right = 0.2
xm = (x_right-x_left)/2.0
f = np.where((x>x_left) & (x<x_right), np.sin(np.pi*(x-x_left)/(x_right-x_left))**2.0,f)
return f
def ftbs_linear(u):
u_new = np.zeros_like(u[1:-1])
u_new[:] = u[1:-1] - (a_max*dt/dx)*(u[1:-1] - u[:-2])
return u_new
def ftbs_nonLin_nonCons(u): # forward time backward space
u_new = np.zeros_like(u[1:-1])
a = wavespeed(u[1:-1]) # wavespeed that varies with u. It is a vector
u_new[:] = u[1:-1] - (a*dt/dx)*(u[1:-1] - u[:-2]) #(1-(dt/dx)*a[1:-1])*u[1:-1] + (dt/dx)*a[1:-1]*u[:-2]
return u_new
def ftbs_nonLin_Cons(u): # forward time backward space
u_new = np.zeros_like(u[1:-1])
Fj = Flux_nonLin(u[1:-1])
Fj_minus = Flux_nonLin(u[:-2])
u_new[:] = u[1:-1] - (dt/dx)*(Fj - Fj_minus) #(1-(dt/dx)*a[1:-1])*u[1:-1] + (dt/dx)*a[1:-1]*u[:-2]
return u_new
def wavespeed(u):
# non-linear wavespeed
a = 0.9 + 0.1*u
return a
def Flux_nonLin(u):
F = 0.9*u + 0.5*0.1*u**2
return F
a_max = 1.0 # linear wave speed
tmin, tmax = 0.0, 1.2 # start and stop time of simulation
xmin, xmax = 0.0, 1.0 # start and end of spatial domain
Nx = 400 # number of spatial points
c = 1.0 # courant number, need c<=1 for stability
# discretize
x = np.linspace(xmin, xmax, Nx+1) # discretization of space
dx = float((xmax-xmin)/Nx) # spatial step size
dt = c/a_max*dx # stable time step calculated from stability requirement
Nt = int((tmax-tmin)/dt) # number of time steps
time = np.linspace(tmin, tmax, Nt) # discretization of time
# initialize variables
u0 = f(x)
u = np.zeros_like(u0)
u[:] = u0[:]
uanalytical = np.zeros((len(time), len(u))) # for storage of the analytical solution
un = np.zeros((len(time), len(u))) # for storage of the numerical solution
# solve from tmin to tmax
#solvers = [ftbs, lax_wendroff, lax_friedrich, macCormack]
solvers = [ftbs_nonLin_nonCons, ftbs_nonLin_Cons]
legends = []
u_solutions = np.zeros((len(solvers), len(time), len(x)))
uanalytical = np.zeros((len(time), len(x))) # holds the analytical solution
for k, solver in enumerate(solvers):
u = f(x)
un = np.zeros((len(time), len(x))) # holds the numerical solution
for i, t in enumerate(time[1:]):
if k == 0:
uanalytical[i + 1,:] = f(x-a_max*t) # compute analytical solution for this time step
u_bc = interpolate.interp1d(x[-2:], u[-2:]) # interplate at right boundary
u[1:-1] = solver(u[:]) # calculate numerical solution of interior
u[-1] = u_bc(x[-1] - wavespeed(u[-1])*dt) # interpolate along a characteristic to find the boundary value
un[i + 1,:] = u[:] # storing the solution for plotting
u_solutions[k,:,:] = un
### Animation
# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure()
ax = plt.axes(xlim=(xmin,xmax), ylim=(np.min(un), np.max(un)*1.1))
lines=[] # list for plot lines for solvers and analytical solutions
legends=[] # list for legends for solvers and analytical solutions
for solver in solvers:
line, = ax.plot([], [], '-', lw=2)
lines.append(line)
legends.append(solver.func_name)
line, = ax.plot([], [], 'k--', lw=2) #add extra plot line for analytical solution
lines.append(line)
legends.append('Analytical (linear advection eq.)')
plt.xlabel('x-coordinate [-]')
plt.ylabel('Amplitude [-]')
plt.legend(legends, loc=3, frameon=False)
# initialization function: plot the background of each frame
def init():
for line in lines:
line.set_data([], [])
return lines,
# animation function. This is called sequentially
def animate(i):
for k, line in enumerate(lines):
if (k==len(solvers)):
line.set_data(x, uanalytical[i,:])
else:
line.set_data(x, u_solutions[k, i,:])
return lines,
# call the animator. blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate, init_func=init, frames=Nt, interval=10, blit=False)
plt.show()