Implementation of 1D Advection in Python using WENO and ENO schemes [closed]

I'm trying to implement 1D advection solver using WENO and ENO schemes.

$$\frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} =0$$

where: $$f(u) = C\,u$$

Discretized as:

$$\frac{u^{n+1}_j-u^{n}_j}{\Delta t} + \frac{F_{j+1/2}-F_{j-1/2}}{\Delta x} =0$$

I've translated Jan S Hesthaven Matlab implementations into Python.

This is how we call the ENO and WENO routines

from numpy import *
nx = 81
dx = 2./(nx-1)
x = linspace(0,2,nx)
nt = 25
dt = .02
c = 1.      #assume wavespeed of c = 1
u = zeros(nx)      #numpy function ones()
u[.5/dx : 1/dx+1]=2  #setting u = 2 between 0.5 and 1 as per our I.C.s
k = 3 # number of weights Order= 2*k-1
gc = k-1 #number of ghost cells
gcr=x[-1]+linspace(1,gc,gc)*dx
gcl=x[0]+linspace(-gc,-1,gc)*dx
xc = append(x,gcr)
xc = append(gcl,xc)
uc = append(u,u[-gc:])
uc = append(u[0:gc],uc)

for n in range(1,nt):
un = uc.copy()
for i in range(1,nx):
xloc = xc[i-(k-1):i+k]
floc = c*uc[i-(k-1):i+k]
#f_left,f_right = ENO(xloc,floc,k)
f_left,f_right = WENO(xloc,floc,k)
uc[i] = un[i]-dt/dx*(f_right-f_left)


Finally either for WENO or ENO I get the same initial condition

This is what we should get using a naive Euler integration.

for n in range(1,nt):
un = u.copy()
for i in range(1,nx):
u[i] = un[i]-c*dt/dx*(un[i]-un[i-1])


Since $\Delta x = (n_t * \Delta t) C = (25*0.02)*1 = 0.5$ the wave displaces 0.5 units to the right.

How should I use the fluxes gotten from WENO and ENO for a simple 1st order Euler integration of 1D advection equation? or Do you have an implementation of 1D Advection using WENO or ENO schemes?

This is the way I'm doing it, but the fluxes ($F_{j+1/2}=$ f_right and $F_{j-1/2}=$f_left) are too small to produce any update. It seems I need to use upwinding for the fluxes.

$$u_{j}^{n+1} = u_j^{n}-\frac{\Delta t}{\Delta x}(F_{j+1/2}-F_{j-1/2})$$

p.s.: For reference only, these are the ENO and WENO routines I'm using:

def ENOweights(k,r):
#Purpose: compute weights c_rk in ENO expansion
# v_[i+1/2] = \sum_[j=0]^[k-1] c_[rj] v_[i-r+j]
#where k = order and r = shift

c = zeros(k)

for j in range(0,k):
de3 = 0.
for m in range(j+1,k+1):
#compute denominator
de2 = 0.
for l in range(0,k+1):
#print 'de2:',de2
if l is not m:
de1 = 1.
for q in range(0,k+1):
#print 'de1:',de1
if (q is not m) and (q is not l):
de1 = de1*(r-q+1)

de2 = de2 + de1

#compute numerator
de1 = 1.
for l in range(0,k+1):
if (l is not m):
de1 = de1*(m-l)

de3 = de3 + de2/de1

c[j] = de3

return c

def nddp(X,Y):
#Newton's divided difference table
#the input are two vectors X and Y that represent points

n = len(X)

DD = zeros((n,n+1))

#inserting x into 1st column of DD-table
DD[:,0]=X

#inserting y into 2nd column of DD-table
DD[:,1]=Y

#creates divided difference coefficients
#e.g: D[0,0] = (Y[1]-Y[0])/(X[1]-X[0])

for j in range(0,n-1):
for k in range(0,n-j-1): #j goes from 0 to n-2
DD[k,j+2]= (DD[k+1,j+1]-DD[k,j+1])/(DD[k+j+1,0]-DD[k,0])

return DD

def ENO(xloc, uloc, k):
#Purpose: compute the left and right cell interface values using an ENO
#Approach based on 2k-1 long vectors uloc with cell k

#treat special case of k=1 - no stencil to select
if (k==1):
ul = uloc[0]
ur = uloc[0]

#Apply ENO procedure
S = zeros(k,dtype=int)
S[0] = k
for kk in range (0,k-1):
#print 'S:',S
#left stencil
xvec = zeros(k)
uvec = zeros(k)
Sindxl = append(S[0]-1, S[0:kk+1])-1
xvec = xloc[Sindxl]
uvec = uloc[Sindxl]
DDl = nddp(xvec,uvec)
Vl = abs(DDl[0,kk+2])

#right stencil
xvec = zeros(k)
uvec = zeros(k)
Sindxr = append(S[0:kk+1], S[kk]+1)-1
xvec = xloc[Sindxr]
uvec = uloc[Sindxr]
DDr = nddp(xvec,uvec)
Vr = abs(DDr[0,kk+2])

#choose stencil through divided differences
if (Vr>Vl):
#print 'Vr>Vl'
S[0:kk+2] = Sindxl+1
else:
S[0:kk+2] = Sindxr+1

#Compute stencil shift 'r'
r = k - S[0]

#Compute weights for stencil
cr = ENOweights(k,r)
cl = ENOweights(k,r-1)

#Compute cell interface values
ur = 0
ul = 0
for i in range(0,k):
ur = ur + cr[i]*uloc[S[i]-1]
ul = ul + cl[i]*uloc[S[i]-1]

return (ul,ur)

def WENO(xloc, uloc, k):
#Purpose: compute the left and right cell interface values using ENO
#approach based on 2k-1 long vectors uloc with cell k

#treat special case of k = 1 no stencil to select
if (k==1):
ul = uloc[0]
ur = uloc[1]

#Apply WENO procedure
alphal = zeros(k)
alphar = zeros(k)
omegal = zeros(k)
omegar = zeros(k)
beta = zeros(k)
d = zeros(k)
vareps= 1e-6

#Compute k values of xl and xr based on different stencils
ulr = zeros(k)
urr = zeros(k)

for r in  range(0,k):
cr = ENOweights(k,r)
cl = ENOweights(k,r-1)

for i in range(0,k):
urr[r] = urr[r] + cr[i]*uloc[k-r+i-1]
ulr[r] = ulr[r] + cl[i]*uloc[k-r+i-1]

#setup WENO coefficients for different orders -2k-1
if (k==2):
d[0]=2/3.
d[1]=1/3.
beta[0] = (uloc[2]-uloc[1])**2
beta[1] = (uloc[1]-uloc[0])**2

if(k==3):
d[0] = 3/10.
d[1] = 3/5.
d[2] = 1/10.
beta[0] = 13/12.*(uloc[2]-2*uloc[3]+uloc[4])**2 + 1/4.*(3*uloc[2]-4*uloc[3]+uloc[4])**2
beta[1] = 13/12.*(uloc[1]-2*uloc[2]+uloc[3])**2 + 1/4.*(uloc[1]-uloc[3])**2
beta[2] = 13/12.*(uloc[0]-2*uloc[1]+uloc[2])**2 + 1/4.*(3*uloc[2]-4*uloc[1]+uloc[0])**2

#compute alpha parameters
for r in range(0,k):
alphar[r] = d[r]/(vareps+beta[r])**2
alphal[r] = d[k-r-1]/(vareps+beta[r])**2

#Compute WENO weights parameters
for r in range(0,k):
omegal[r] = alphal[r]/alphal.sum()
omegar[r] = alphar[r]/alphar.sum()

#Compute cell interface values
ul = 0
ur = 0
for r in range(0,k):
ul = ul + omegal[r]*ulr[r]
ur = ur + omegar[r]*urr[r]

return (ul,ur)

• Presumably you've tried to solve any issues on your own. At a minimum, please let us know what you've tried and any clues you have as to what's wrong. Does the code run or are you getting errors? Are the results wrong? If so, in what way? – Doug Lipinski Jul 1 '15 at 17:29
• Yes I've compared the Python WENO and ENO subroutines with matlab, so the routines should be OK. I'm not sure if the euler integration is correct. The values for the fluxes are too small to produce any advection. – ilciavo Jul 1 '15 at 17:52
• You mean Jan Hesthaven (not "Jed"). – David Ketcheson Jul 2 '15 at 8:31
• I downvoted because you're simply asking others to do the work of debugging your code. – David Ketcheson Jul 2 '15 at 8:32
• It's not really clear what you're asking. You've posted some code and then said "Let me know if you can spot my mistake." To me, this feels like a scavenger hunt without clues. You've implied that there is a mistake, but not given us any idea how you know there is a mistake, what the mistake might be, or even if the mistake is in your code or in your concepts. Like @DavidKetcheson, I also thought you were just asking for code debugging. Please clarify your question. – Doug Lipinski Jul 2 '15 at 12:49

The confusion was the misleading variables $F_{j-1/2}$ and $F_{j+1/2}$ with f_left and f_right, which are completely different.

f_left and f_right are the interpolated fluxes at a one single face. They must be then upwinded using the advection speed to compute the Flux at a specific cell face. Which means if $C>0$ we take f_left, otherwise we take f_right.

gs = zeros((nx+2*gc,nt))
flux = zeros(nx+2*gc)

for n in range(1,nt):
un = uc.copy()
for i in range(gc,nx-1+gc): #i=2
xloc = xc[i-(k-1):i+k] #i+k-1-(i-(k-1)-1) = 2k -1
uloc = uc[i-(k-1):i+k]
f_left,f_right = ENO(xloc,uloc,k)
#f_left,f_right = WENO(xloc,uloc,k)
#upwind flux
flux[i]=0.5*(c+fabs(c))*f_left + 0.5*(c-fabs(c))*f_right

for i in range(gc,nx-gc):
if c>0:
uc[i] = un[i]-dt/dx*(flux[i]-flux[i-1])
else:
uc[i] = un[i]-dt/dx*(flux[i+1]-flux[i])


Finally this is the advected scalar