# What is wrong in the code for this upwind method?

I want to implement the upwind method in following advection equation problem :
$$u_{t}+2 u_{x} =0 ,$$ for $$0\leq t \leq 1,$$ $$0\leq x \leq 1$$ $$u(0,x) = u_{0}(x) = \begin{cases} 10^4 (0.1-x)^2 (0.2-x)^2 & 0 < x < 0.2 \\ 0 & otherwise \end{cases}$$ with wind coefficient $$\alpha = 2$$ , $$Nt = 100$$ ,$$Nx = 250$$ resulting to Courant number $$v = 0.4$$

But it doesn't seem right.What i do wrong? I would appreciate any help.

import numpy as np
import matplotlib.pyplot as plt

class UpwindMethod:

def __init__(self, Nx,Nt, tmax):
self.Nx = Nx                # number of space nodes
self.Nt = Nt                # number of time  nodes
self.tmax = tmax
self.xmin = 0
self.xmax = 1
self.dt =   tmax/self.Nt    # timestep
self.a =    2               # velocity
self.v =    (self.a * self.dt ) / ((self.xmax - self.xmin)/self.Nx)
self.initializeDomain()
self.initializeU()
self.initializeParams()

def initializeDomain(self):
self.dx = (self.xmax - self.xmin)/self.Nx
self.x = np.arange(self.xmin-self.dx, self.xmax+(2*self.dx), self.dx)

def initializeU(self):
u0 = np.where( (0.1<self.x) & (self.x<0.2) , ((10**4)*((0.1-self.x))**2 *((0.2-self.x)**2))   ,0)
self.u = u0.copy()
self.unp1 = u0.copy()

def initializeParams(self):
self.nsteps = round(self.tmax/self.dt)

def solve_and_plot(self):
tc = 0.

for n in range(self.nsteps):
plt.clf()
for i in range(self.Nx+2):
self.unp1[i] = self.u[i] - ((self.a*self.dt)/(self.dx))*(self.u[i]-self.u[i-1])

self.u = self.unp1.copy()

# Periodic boundary conditions
self.u[0] = self.u[self.Nx+1]
self.u[self.Nx+2] = self.u[1]

uexact = np.where( (0.1<self.x) & (self.x<0.2) , ((10**4)*((0.1-self.x-self.a* tc)**2) *((0.2-self.x-self.a* tc)**2))   ,0)

plt.plot(self.x, uexact,  label=" Exact solution ")
plt.plot(self.x, self.u,  label=" Upwind Approximation ")

plt.xlabel("x")
plt.ylabel("u")
plt.legend(loc=1, fontsize=13)
tc += self.dt

def main():
sim = UpwindMethod(100,250, 0.02)
sim.solve_and_plot()
plt.show()
if __name__ == "__main__":
main()

$$$$

• It might help to show a plot of the solution and discuss why you think it is incorrect. – Paul Nov 13 '20 at 14:49

I've identified a few "problems":

1. Your analytical solution isn't quite correct. The correct analytical solution to the "infinite domain" advection equation is supposed to be $$u(t,x) = u_0(x-at)$$. Because you have periodic boundaries, you need to properly account for this by making sure $$x-at$$ is properly wrapped.

I replaced your initializeU function with one which just computes u(t,x):

def exact_u(self, t):
# this properly wraps your analytical solution for periodic BC's
xp = (self.x - self.a*t - self.xmin)%(self.xmax-self.xmin) + self.xmin
u0 = np.where( (0.1<xp) & (xp<0.2) , ((10**4)*((0.1-xp))**2 *((0.2-xp)**2)), 0)
return u0


You will also need to update the __init__ or initializeU function to correctly use this function by just calling it for the exact solution at $$t=0$$.

In addition to this, where you're using this in the solve_and_plot function is comparing the simulation results with the exact solution 1 timestep prior, so you should change it to:

tc += self.dt
uexact = self.exact_u(tc)

1. I don't think you're handling the periodic BC's correctly in your simulation. There's not really any need to create extra nodes at the exterior, since you can take advantage of Python's negative indexing to index to the end of the array automatically with your 2 node upwind stencil. You also have to take into account the fact that the very first node and last node are effectively "coincident" and supposed to be identical.

You can handle this by stopping one $$\Delta x$$ short at the end, since it's just a duplicate node, however this does look a bit weird since you're not plotting to the end of your periodic domain. Here's a solution which circumvents this by doing some slicing to handle the simulation correctly while still having an extra storage space for the end node:

    def initializeDomain(self):
self.dx = (self.xmax - self.xmin)/self.Nx
self.x = np.linspace(self.xmin, self.xmax, self.Nx + 1)

def solve_and_plot(self):
tc = 0.

for n in range(self.nsteps):
plt.clf()
# periodic BC's automatically handled by Python's negative indexing
# use slices to ensure we skip the duplicate end node only there for plotting purposes
unp1 = self.unp1[:-1]
u = self.u[:-1]
for i in range(u.shape[0]):
unp1[i] = u[i] - ((self.a*self.dt)/(self.dx))*(u[i]-u[i-1])

# update the duplicate end node's value
u[-1] = u[0]

self.u = self.unp1.copy()

#... rest of this function

1. As you've probably noticed when you have a $$v < 1$$, the simulation is stable but diffusive (becoming less diffusive as $$v \rightarrow 0$$. This is simply a fundamental numerical error associated with how you chose to discretize the method. You can "magically" avoid this for this problem by using a timestep which gives you exactly a $$v = 1$$ (this doesn't always work). In this case, this requires $$\Delta t = \frac{\Delta x}{a}$$, or Nt = 4`. As a side-note, I think you're calculating your Courant number incorrectly; it's supposed to be: $$v = \frac{a \Delta t}{\Delta x}$$

$$v = 1$$ (Nt = 4, $$t_{end} = 0.02$$):

$$v = 0.016$$ (Nt = 250, $$t_{end} = 0.02$$):