How to implement SUPG properly for advection dominated equation

Question

Hi All,

I am trying to solve a transient advection problem, and would like to know how to properly implement the SUPG method to stabilize my solution. Thus, I use the SUPG example from fenics as my starting point, and try to modify it.

It is to my understanding that different people have proposed different formulations for $\tau$, which is a time scale scalar. In SUPG fenics example, tau = h/(2.0vnorm), and I have read a few other textbooks indicating more complicated terms in addition to the h/(2.0vnorm). (equation added in the end)

Here is my test case to see if I understand it correctly, I have two tests,

1. Applied a Dirichlet sine function boundary, and let this sine function move with respect to the velocity field.
2. Initiate a concentration of 1 on a circle, and let it move with velocity field, and see how SUPG can keep the shape of the circle.

For the test 1, both tau appears to behavior ok. However, it surprise me that the more complicated tau that from the textbook always have higher overshoot (highest goes to 1.15, and lowest goes down to -1.15), while the more simple one, is always within 1.1 and -1.1.

For the test 2, SUPG based on these two tau behave pretty bad, and introduce further instability as the velocity moves the circle.

My question is: First, did I implement SUPG properly ( I have attached my full code in below) ? Second, If I implement properly, then is this what should be expected from SUPG method ? From what I read, SUPG should perform much better than what I see. Third, did I miss any important stabilization criterion that could lead to the stabilization I see in test 2?

Thanks

from dolfin import *

# Create mesh and function space
mesh = UnitSquareMesh(100, 100)
Q = FunctionSpace(mesh, "CG", 1)

# Define velocity as vector
velocity = Expression(("1" , "1"), domain = mesh, degree=1)

c = 1e-15
dt = 0.01
# Test and trial functions
u, v = TrialFunction(Q), TestFunction(Q)
u0 = Function(Q)

#initial condition
ic= Expression("((pow(x[0]-0.2,2)+pow(x[1]-0.2,2))<0.1*0.1)?(1.0):(0.0)", domain=mesh, degree = 3)
u0.interpolate(ic)

# Crank Nicolson
u_mid = 0.5*(u0 + u)

# Residual
r = u - u0 + dt*(dot(velocity, grad(u_mid)) - c*div(grad(u_mid)))

# Galerkin variational problem
F = v*(u - u0)*dx + dt*(v*dot(velocity, grad(u_mid))*dx + c*dot(grad(v), grad(u_mid))*dx)

# Add SUPG stabilisation terms
vnorm = sqrt(dot(velocity, velocity))
h = CellSize(mesh)
tau = h/(2.0*vnorm) # tau from SUPG fenics example
#tau = pow(1/(0.5*dt) + 2.0*vnorm/h + 4*c/pow(h,2.0),-1) # tau from
F += tau*dot(velocity, grad(v))*r*dx

# Create bilinear and linear forms
a = lhs(F)
L = rhs(F)

# Define boundary condition
tol = 1e-14
def u0_boundary(x, on_boundary):
return on_boundary and near(x[1], 0, tol)
U0_expression = Expression('sin(x[0]*60)', domain = mesh, degree=3)
#bc = DirichletBC(Q, U0_expression, u0_boundary)
vtkfile = File('results/concentration.pvd')

t = 0;
T = 1.0;
while t < T:
# Solve
    u = Function(Q)
    #solve(a == L, u, bc)
    solve(a == L, u)
    u.rename("u", "concentration")
    vtkfile << u
    u0.assign(u)
    t = t + dt

Comments

How does the textbook define h? If I'm not mistaken FEniCS used the cell diameter, but that is one of many suitable definitions and the cell shape and orientation w.r.t. the flow may play a role here.

---

Well, compared to the pure Galerkin formulation I guess SUPG performs not to badly in case 2? Did you take the problem setup from the literature, is there a reference solution? Otherwise I'd suggest to do a well studied and documented setup first and compare your results with the literature.

---

Hi, the textbood define h as meshSize, which is not very helpful on whether it is radius or diameter. However, in the introduction of SUPG, the textbook introduced a critical tau for 1D convection-diffusion as tau = h/2*a (coth(Pe) - 1/Pe), which is exactly the same as used in fenics 2011 book (Page 637). So I assume the h from the textbook is consistent with what is implemented in fenics.

---

Hi dajuno, Yes. I agree that SUPG performs better than the pure Gakerkin. I probably should have done more detailed literature setup study. However, It just sometimes surprise me that a more completed terms which its self contains the simple version of the tau and several other many other terms perform "worse" than the simply tau = h/2a. My initial thought about this is that the the simple tau expression will basically be the base line of SUPG, and more complicated tau would do better.

That is why I wonder whether myself implement this correctly or not....

---

Well, I guess it just depends. The implementation of SUPG looks ok, don't know if the test case is "good".
Did you compare the values of the $\tau$s? It seems that in your setup ($c\to 0$, $\Delta t$ "large") the formulas are pretty much identical.

---

For problems like case 2, with (nearly) discontinuous solutions, you might look into discontinuity capturing methods, e.g.

https://www.ices.utexas.edu/media/reports/2007/0702.pdf

---

Hi Kamensky, sorry for the late reply. Thanks for pointing out additional literature resources. I enjoyed reading it. It is interesting to see the approach using artificial viscosity on top of SUPG. I will try to implement it. Just curious, do you have any working fenics example similar to the approach described in the paper that is available to share ?