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

Question

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)

Answer 1

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