I'm trying to implement 1D advection solver using WENO and ENO schemes.
\begin{equation} \frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} =0 \end{equation}
where: \begin{equation} f(u) = C\,u \end{equation}
Discretized as:
\begin{equation} \frac{u^{n+1}_j-u^{n}_j}{\Delta t} + \frac{F_{j+1/2}-F_{j-1/2}}{\Delta x} =0 \end{equation}
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
#adding 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)