Lax-Friedrichs Flux for advection equation

1

Hi everyone,

I am currently trying to replace the upwinding flux for the advection equation \dot{u} + \nabla \cdot (\beta u) = f.

       un = (dot(beta, n) + abs(dot(beta,n)))/2.0
       flux = jump(un*phi)

with a Lax-Friedrichs flux:

        C = beta[0]
        flux = avg(dot(beta*phi,n)) + 0.5*C*jump(phi)

I used the definition from Symbolic Computation for DGFEMs where
\mathcal{H}_{LF}(u^+_h, u_h^-, n_k) = \frac{1}{2} (\mathcal{F}(u_h^+)\cdot n_k + \mathcal{F}(u_h^-)\cdot n_k + \alpha (u^+_h - u^-_h) ) \quad (3.1a).
For this simple test case the maximum eigenvalue \alpha should simply be the wave speed in x-direction.

However the solution computed using the this definition is clearly wrong as the it becomes zero.

Does anyone perhaps see what I am formulating incorrectly in my code? Thank you in advance! :slightly_smiling_face: As the definition is really quite straight forward, I just don’t see my error at the moment.

My complete code is:

def forward(beta, mesh, delta_t = 0.01, end_time=1, flux_method = 'upwind', show_plots=False):
    parameters["ghost_mode"] = "shared_facet"
    
    V_dg = FunctionSpace(mesh, 'DG', 2)
    phi, v = TrialFunction(V_dg), TestFunction(V_dg)
    u_new, u_old = Function(V_dg), Function(V_dg)
    
    u0 = Expression("sin(2*pi*x[0])", degree=1, t=0, name='u0') # 
    u_old = project(u0, V_dg)
    u_exact = Expression("sin(2*pi*(x[0] - c*t))", degree=2, t=0, c=beta[0]) # 
    bcs = DirichletBC(V_dg, u_exact, "on_boundary", method='geometric')
    
    f = Constant(0.0) 
    time_steps = int(end_time/delta_t)
    
    n = FacetNormal(mesh)
    F_u = beta*phi

    if 'LaxF' == flux_method:
        C = beta[0]
        flux = avg(dot(beta*phi,n)) + 0.5*C*jump(phi)
    elif 'upwind' == flux_method:
        un = (dot(beta, n) + abs(dot(beta,n)))/2.0
        flux = jump(un*phi)

    F = Constant(1/delta_t)*(phi-u_old)*v*dx - dot(grad(v), F_u)*dx -  v*f*dx
    F += dot(jump(v), flux)*dS + dot(v, u_exact*phi)*ds
    
    a,L = lhs(F),rhs(F)

    t = 0
    for n in range(time_steps):        
        solve(a==L, u_new, bcs)
        u_old.assign(u_new)
        
        t += delta_t
        u_exact.t = t
        
    if show_plots:
        plot(u_old)
        plt.show()
        
    return u_old

# compare upwind flux with lax-friedrichs flux
m_size = 32
mesh = UnitSquareMesh(m_size, m_size)
beta_x = Constant(0.5) 
beta_y = Constant(0) 
beta = as_vector((beta_x, beta_y))

u_up = forward(beta, mesh, flux_method='upwind')
u_lf = forward(beta, mesh, flux_method='LaxF')

plot(u_up)
plt.show()

plot(u_lf)
plt.show()
1 Like
2

I don’t really have time to look at your code in detail, but your implementation of the local-LF flux is a little confusing. One verification tip would be to ensure your symbolic formulation for the local-LF flux and the upwind flux is equivalent, as it should be in the linear setting.

1 Like
3

It looks to me like a misunderstanding of how normals work on interior facets. If I make the replacement

        #flux = avg(dot(beta*phi,n)) + 0.5*C*jump(phi)
        flux = dot(avg(beta*phi),n('+')) + 0.5*C*jump(phi)

the L–F flux produces similar results to upwinding. This is consistent with the interpretation of \mathbf{n}_\kappa as n('+') from Listing 2 from @nate’s paper.

3 Likes
4

Thank you for looking so closly @kamensky! Such an obvious mistake on my end.

@nate Thanks for the debugging tip!

1 Like
6

Hi, thank you for this discusson. I can’t help but wonder why there is no n(‘+’) at the end of 0.5Cjump(phi).

The jump of phi (which is a scalar) would be a vector and shouldn’t we be having its dot product with the normal?