Neumann boundary conditions for the upwind scheme applied to the advection equation (Python)

Question

I'm trying to solve the linear advection equation

$$u_{t} = cu_{x}, \\ x \in [x_{0}, x_{e}], \quad t \in (0, T], \quad c \in \mathbb{R} \\ u(x,0) = f(x)$$

Note that for $c > 0$, the solution is a left travelling wave so requires a boundary condition on the left, while $c < 0$ is a right travelling wave so requires a boundary condition on the right. Now, implementing Dirichlet ($u(\cdot, t) = 0$) and periodic ($u(x_{0},t) = u(x_{e},t)$) boundary conditions has been quite straight forward. However, I've been trying to implement Neumann conditions ($u_{x}(\cdot,t) = \alpha$) and I'm not exactly sure how to or if it is possible?

If $c > 0$, then the discretisation in matrix form is given as

$$\begin{pmatrix} u_{1}^{k+1} \\ u_{2}^{k+1} \\ u_{3}^{k+1} \\ \vdots \\ u_{N}^{k+1} \end{pmatrix} = \begin{pmatrix} 1+\gamma \\ -\gamma & 1+\gamma \\ & -\gamma & 1+\gamma \\ & & \ddots & \ddots \\ & & & -\gamma & 1+\gamma \end{pmatrix} \begin{pmatrix} u_{1}^{k} \\ u_{2}^{k} \\ u_{3}^{k} \\ \vdots \\ u_{N}^{k} \end{pmatrix}$$

If we are to use the ghost point approach, then applying a central difference we get

\begin{align} u_{x}(x_{N},t) &= \frac{u_{N+1} - u_{N-1}}{2 \Delta x} \\ &= \alpha \\ \implies u_{N+1} &= u_{N-1} + 2 \alpha \Delta x, \quad \forall k \end{align}

The problem is, nowhere in our spatial discretisation up to $n = N$ do we have our function indexed at $N+1$, so we can't use that approach. If we use a backward difference instead, we get

\begin{align} u_{x}(x_{N},t) &= \frac{u_{N} - u_{N-1}}{\Delta x} \\ &= \alpha \\ \implies u_{N} &= u_{N-1} + \alpha \Delta x, \quad \forall k \end{align}

and hence our final ODE becomes

$$u_{N}^{k+1} = u_{N}^{k} + \alpha \Delta x$$

but that didn't seem to work either. So my question is, how do we implement Neumann boundary conditions using the upwind method?

Here is my code

def Upwind(N,x0,xe,t0,te,m,c,alpha,f,boundary='Dirichlet'):

# ----------------------
### Timestep

t = np.linspace(t0, te, N)
dt = t[1]-t[0]

# ----------------------
### Direction conditions

if c == 0:
    
    print('u(x,t) = f(x)')

else:
    
    if c > 0:

        x = np.linspace(xe, x0, N)
        
    else:

        x = np.linspace(x0, xe, N)
        
        
    dx = x[1]-x[0]
    r = c*np.ones(N)*dt/dx

    # ----------------------
    ### Stability conditions

    if abs(r[0]) > 1+1e-15:

        print('r = {}'.format(r[0]))
        print('This solution is unstable, abs(r) must be less \
              than or equal to one for stability.')
           
    # ----------------------
    ### Derivative matrix
    
    DL = spdiags([np.ones(N)], (0), N, N).tocsr() # NOT NECESSARY, MORE FOR 
                                                  # CONTINUITY WITH OTHER SCHEMES 
                                                  # THAT HAVE A LHS
    DR = spdiags([-r, 1+r], (-1, 0), N, N).tocsr()

    # ----------------------
    ### Boundary conditions

    if boundary == 'Dirichlet': # AGAIN, NOT NECESSARY, MORE FOR 
                                # CONTINUITY. HERE u(.,t) = 0.

        BL = spdiags([0], (0), N, N)
        BR = spdiags([0], (0), N, N)

    elif boundary == 'Neumann': # NEEDS TO BE FIXED

        BL = spdiags([0], (0), N, N)
        BR = spdiags([0], (0), N, N)

    else:
        
        BL = spdiags([0], (0), N, N)
        BR = spdiags([-r], (N-1), N, N)

    # ----------------------    
    ### Initial condition

    u = f(x)

    # ----------------------
    ### Loop

    for k in range(m+1):

        u = spsolve(DL+BL, (DR+BR)*u)
        
        if k % 25 == 0:
        
            plt.plot(x,u)

and to test

Upwind(200,-40,40,0,40,100,-2,1,lambda x: np.exp(-x**2),boundary='Periodic')

Answer 1

The problem is, nowhere in our spatial discretisation up to n=N do we have our function indexed at N+1

You just add an extra element (i.e. at N+1).

Then, every loop set it to the value of the last element times $\alpha\Delta x$. When the forward difference is calculated, it uses a value at $i+1$ to work out $i$. So although there's an extra element, it isn't part of the domain and doesn't appear in the result: it is purely a condition on the boundary.

Alternatively, in your code, you can handle the forward difference for the last element as a special case: instead of looking up the next element (i.e. past the last element), use the last element times $\alpha\Delta x$.

BTW: I think Barba introduces boundary conditions and other issues more simply in her 12 Steps to Navier Stokes (which don't use a matrix solver). Might be worth a look; though it is more work.

Answer 2

The time-derivative scheme you are using is referred to as Forward Euler as is one of the simplest schemes out there, which I will assume for the rest of my answer. Even in the case that you would like to extend this, it is quite straightforward.

Coming to your implementation, we basically need to use \begin{equation} u_{j}^{N+1} = u_{j}^{N} + c \Delta t u_{x} \end{equation}

where you replace $u_{x}$ with the numerical derivative at the interior points and with the boundary condition for the end points.

Given that you've stuck to a matrix notation, which complicates the implementation, the only change you'll need to make is

RHS = (DR + BR)*u
if boundary == 'Neumann':
    RHS[0] += c*dt*(alpha-u[0]/dx)
u = spsolve(DL+BL, RHS)

This replaces the left boundary derivative with the boundary condition.

The Wikipedia article for the first-order upwind scheme is clear and would recommend going through it (https://en.wikipedia.org/wiki/Upwind_scheme).

Answer 3

You seem to be using explicit Euler in time, so your scheme would be \begin{equation} u^{k+1}_j = u^k_j + \frac{c \Delta t}{\Delta x} (u^k_j-u^k_{j-1}) \end{equation} (if so, I don't really understand why you would make a sparse matrix). With the second approach that would give \begin{equation} u^{k+1}_N = u^k_N + \alpha c \Delta t \end{equation}