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 
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.

Naive Euler

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 

    #inserting y into 2nd column of DD-table

    #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
            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):
        beta[0] = (uloc[2]-uloc[1])**2
        beta[1] = (uloc[1]-uloc[0])**2

        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)
  • $\begingroup$ 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? $\endgroup$ – Doug Lipinski Jul 1 '15 at 17:29
  • $\begingroup$ 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. $\endgroup$ – ilciavo Jul 1 '15 at 17:52
  • $\begingroup$ You mean Jan Hesthaven (not "Jed"). $\endgroup$ – David Ketcheson Jul 2 '15 at 8:31
  • $\begingroup$ I downvoted because you're simply asking others to do the work of debugging your code. $\endgroup$ – David Ketcheson Jul 2 '15 at 8:32
  • 1
    $\begingroup$ 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. $\endgroup$ – 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])
            uc[i] = un[i]-dt/dx*(flux[i+1]-flux[i])

Finally this is the advected scalar

enter image description here

| cite | improve this answer | |

Not the answer you're looking for? Browse other questions tagged or ask your own question.