Burger's equation (PDE) does not work with downwind difference?

Question

I'm working on implementing the discretised Burger's equation. I am quite confused as to why it does not work when using a step-function and downwind difference formula. When using a step-function and downwind difference, Newton's method fails to converge.

Utilising downwind difference should provide results, as using upwind seemed to work. but I'm unable to identify where or why it's going wrong. I'm trying to run it with:

viscosity value ν = 0.001

time-step t = 0.01

grid size h = 0.01

and final time T = 1.

The initial solution u(x, 0) is taken as a unit step at x = 0.1, boundary conditions as: u(0, t) = 1 and u(1, t) = 0.

Would greatly appreciate insights into where I am going wrong.

I've provided a reproducible example: burgers_solve.py

This is how I'm implementing the discretised Burger's equation with downwind:

# Downwind
def fun_burgers(uk, ukp, dt, h, nu, ua, ub):
    # ukp is previous time step solution
    # dt is time-step t
    # h is spatial grid size parameter 
    # nu is kinematic viscosity 

m = ukp.size
# f to store values of the function for each point uk in space
f = np.zeros((m-2, 1))

# boundary conditions
uL = ua
uR = ub

# Left boubdary
f[0] = (uk[0] - ukp[1])/dt + uk[0] * (uk[1] - uk[0])/h - nu * (uk[1]-2*uk[0] + uL)/h**2

# Difference equations at each internal node
for i in range(1, m-3):
    f[i] = (uk[i] - ukp[i+1])/dt + uk[i] * (uk[i+1] - uk[i])/h - nu * (uk[i+1] - 2*uk[i] + uk[i-1])/h**2

# Right boundary
f[m-3] = (uk[m-3] - ukp[m-2])/dt + uk[m-3] * (uR - uk[m-3])/h - nu * (uR - 2*uk[m-3] + uk[m-4])/h**2

return f

Answer 1

The qualitative difference between upwind and downwind flux (note that both are first order accurate for twice continuously differentiable solutions) lies in the fact that the upwind flux (in the context of finite volume schemes)

\begin{align} F_{i + \frac{1}{2}} = \begin{cases} \min_{U_i^{(n)} \leq \theta \leq U_{i+1}^{(n)}} f(\theta) & U_i^{(n)} \leq U_{i+1}^{(n)} \\ \min_{U_{i+1}^{(n)} \leq \theta \leq U_{i}^{(n)} } f(\theta) & U_i^{(n)} > U_{i+1}^{(n)} \end{cases} \end{align}

actually corresponds to the true solution of the Riemann problem (see for instance Chapter 2.1 of Godlewski & Raviart) \begin{align} \partial_t u + \partial_x u^2 &=0 \\ u\Big(x, t^{(n)} \Big) &= \begin{cases} u_i^{(n)} & x < x_{i + \frac{1}{2}} \\ u_{i+1}^{(n)} & x > x_{i + \frac{1}{2}} \end{cases} \end{align}

at the cell interfaces $x_{i \pm \frac12}$.

You can even show for the simpler advection equation $$u_t + a u_x = 0, \quad a = \text{const} $$ with periodic boundary conditions that the usage of the downwind discretization leads to an unbounded growth of the (numerical) energy of the solution

$$E^{(n)} = \frac{1}{2} h \sum_i \left( u_i^{(n)} \right)^2 $$ which indicates more trouble for non-constant coefficient systems such as Burgers' equation.