Convergence Issues for 1D Burgers Equation

variational formulation
1

To preface, I’m very new to Fenics (& finite elements) so I appreciate any help or pointers I can get. I’m working on solving Burgers Equation, but I’m running into some issues achieving convergence.

This is the problem I’m working on:
u_t + uu_x - (0.01/\pi)u_{xx}= 0
u(0, x)=-sin(\pi x)
u(t,-1)=u(t,1)=0
where x \in [-1,1] and t \in [0,1].

I’ve defined the variational formulation as follows:
F = [\frac{u-u_n}{dt}\cdot v + u\nabla u \cdot v - \frac{1}{100\pi}\nabla u\nabla v]dx

I’ve messed with the parameters for a while now, but is there anyway I could better define the variational formulation? I’ve messed around with the number of mesh points and the number of time steps without too much luck. I’ve also tried using a few different solvers instead of the default NewtonSolver as well, but not much luck there either.

Here’s a working example of what I have so far:

from fenics import *
import matplotlib.pyplot as plt
%matplotlib notebook

nx = 100
x0, xf = (-1.0, 1.0)
T = 1.0 
nSteps = 1000
dt = T/nSteps
order = 2

mesh = IntervalMesh(nx, x0, xf)
V = FunctionSpace(mesh, "CG", order)

# boundary conditions
def left(x, on_boundary):
    return on_boundary and near(x[0], x0)
def right(x, on_boundary):
    return on_boundary and near(x[0], xf)

bc_l = DirichletBC(V, Constant(0), left)
bc_r = DirichletBC(V, Constant(0), right)
bcs = [bc_l, bc_r]

# initial conditions
u_0 = Expression('-sin(pi*x[0])', degree=1)
u_n = interpolate(u_0, V)

# define variational problem
u = Function(V)
v = TestFunction(V)
u_x = u.dx(0)

DT = Constant(dt)
alpha = Constant(1/(100*pi))
f = Expression("0.0", degree=0)

F = (dot(u - u_n, v)/DT
    + dot(u*u_x, v)
    - alpha*dot(grad(u), grad(v))
    - dot(f, v))*dx

t = 0
for n in range(nSteps):
    
    # update time
    t += dt

    # solve variational problem
    solve(F == 0, u, bcs)
    
    # update previous solution
    u_n.assign(u)

This seems like a pretty simple problem to solve, so I think I must be setting it up incorrectly somehow. Thanks for the help.

2

The immediate problem is that you forgot to flip the sign when integrating by parts in the viscous term, so your weak formulation corresponds to unstable anti-diffusion.

Fixing that, the nonlinear convergence is fine, but I’ll mention that the Bubnov–Galerkin method is not robust in the inviscid limit of \alpha\to 0. (In general, you’d need \vert u\vert h/\alpha \lesssim 1 for good results, where h is the element size, implying that a very fine mesh is needed for high Reynolds numbers.) You can improve on this using SUPG stabilization, e.g.,

# Advective flux of conservation law:
e0 = Constant((1,))
def flux(u):
    # Vector of dimension 1, for consistency w/ grad, div
    return 0.5*u*u*e0

# Advection velocity for a general hyperbolic conservation law:
u_var = variable(u)
a = diff(flux(u_var),u_var)

# Limit explosion of automatically-determined quadrature degree; just try to get the
# advection term exactly:
dx = dx(metadata={"quadrature_degree":3*order-1})

# Galerkin part of weak problem:
F_Gal = (dot(u - u_n, v)/DT
         - dot(flux(u),grad(v))
         + alpha*dot(grad(u), grad(v))
         - dot(f, v))*dx

# Residual of strong problem:
res_strong = (u-u_n)/DT + div(flux(u) - alpha*grad(u)) - f

# SUPG stabilization:
h = CellDiameter(mesh)
C_inv = Constant(6)
C_t = Constant(2)
tau_adv = h/(2*sqrt(dot(a,a) + DOLFIN_EPS))
tau_diff = h*h/(C_inv*alpha + DOLFIN_EPS)
tau_t = DT/C_t
tau = 1.0/sqrt(tau_adv**(-2) + tau_diff**(-2) + tau_t**(-2))
F_SUPG = tau*res_strong*dot(a,grad(v))*dx

# Combined formulation; comment out `+ F_SUPG` with
# `alpha = Constant(DOLFIN_EPS)` to see effect of stabilization.
F = F_Gal + F_SUPG

#F = (dot(u - u_n, v)/DT
#    + dot(u*u_x, v)
#    + alpha*dot(grad(u), grad(v))
#    - dot(f, v))*dx

“Industrial grade” solvers for hyperbolic conservation laws also usually include some mechanism to add further dissipation around shocks, but shock-capturing viscosities that preserve high-order accuracy for smooth solutions are usually very nonlinear and present major convergence difficulties when used in conjunction with implicit time integrators like this.

For the theory behind SUPG, some good references are the original paper, and a later numerical analysis of the method for advection–diffusion. It’s also covered (again for advection–diffusion) in the lecture notes for my graduate class at UCSD (from the more recent variational multiscale (VMS) perspective). There are probably some good references specifically on SUPG for the Burgers equation, but I’ve never looked into it before.

4 Likes
3

Wow, duh. I appreciate the help and the additional knowledge about improving stability. Luckily, I won’t need to worry too much about this, as 100\pi is the only Reynold’s number I’ll be working with for my purposes. But nonetheless, it’s still very good to know.

I’m gonna have to do some more digging to fully appreciate your SUPG implementation. I’ll definitely be looking into the sources you mentioned and your course github repo as well.
Thanks!