Basic 1D advection equation in python using the finite difference method

Question

Hi I am trying to code a simple advection equation in python using the finite difference upwind method. I have looked online to find a simple example of this but the codes I have found are a little more involved than what I expected to find. I have a simple code written in python where I have used a Gaussian initial condition. The equation I am working with and code I have written are below:

Equation

import numpy as np
import matplotlib.pyplot as plt

#Mesh
Nx = 211 #Number of x grid points
Nt = 1000 #Number of t grid points

dx = 1/(Nx - 1) #Dividing dx into equal dx intervales
dt = 1e-7 #Time step

#Boundary conditions
x= np.linspace(0.1, 1.0, num=Nx)  
t= np.linspace(0.0, Nt*dt, num=Nt)

#Initializing array to zero
u=np.zeros((Nx, Nt)) #Initializing array to zero

#Velocity
c = 1 

#Gaussian initial condition
mu = 0.5
sigma = 0.1

for i in range(1,Nx-1):
    u[i,0] = np.exp(-(x[i]-mu)**2/(2*sigma**2)) / np.sqrt(2*np.pi*sigma**2)


#Looping over all grid points excluding those at the boundary conditions
for j in range(0,Nt-1):   
    for i in range(1,Nx-1):
        
        u[i+1,j] = u[i,j] + c*(((1/x[i]) + (u[i+1,j] - (u[i,j])/dx)))*dt

plt.plot(x, u[:, 0])
plt.plot(x, u[:, 200])
plt.plot(x, u[:, 300])
plt.plot(x, u[:, 400])

I was expecting the plot to move to the left since the speed c>0 and wanted to see the plot at 4 different times beginning with the initial Gaussian starting point. But what I get doesn't make any sense. I also want to point out that when I comment out the for loop and do plt.plot(x,u[:,0] I see the initial Gaussian curve as expected when I type it again after uncommenting the for loop I get the plot shown below. I would greatly appreciate any help anyone can give me with this.

Output Plot

Answer 1

You can try this. Note, in particular:

I have temporarily switched off the source term (c/r) by multiplying it by 0. An accurate solution (increase Nx) will then advect unchanged.

I have corrected your du/dx term and the time level on the LHS.

With your Forward-Euler scheme you only need TWO time levels, plotting as you go along.

I have increased the number of x-nodes for greater accuracy. I have also made it consistent with dx.

You still need to think about your boundary conditions.

import numpy as np
import matplotlib.pyplot as plt

# Time dependence
Nt = 1000          # number of timesteps
dt = 1e-3

# Space dependence
Nx = 901           # number of nodes
xmin, xmax = 0.1, 1.0
dx = ( xmax - xmin ) / (Nx - 1)
x= np.linspace(xmin, xmax, Nx)

# Equation variable (speed)
c = 1 

# Dependent variable with gaussian initial condition
mu = 0.5
sigma = 0.1
u = np.exp(-(x-mu)**2/(2*sigma**2)) / np.sqrt(2*np.pi*sigma**2)
plt.plot(x, u)

# Loop over all grid points excluding those at the boundary conditions (this IS a problem)
for j in range( 1, Nt + 1 ):
    uold = u.copy()
    for i in range(1, Nx - 1):
        u[i] = uold[i] + c * dt * ( uold[i+1] - uold[i] ) / dx    + 0 * c / x[i]
    if j % 100 == 0: plt.plot( x, u )

plt.show()

Output: