This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Help on SUPG method in the advection-diffusion demo

+2 votes

Hi there,

I'm studying the SUPG method in the advection-diffusion demo, and tested a manufactured stationary problem on a unit square mesh.
Basically, I kept the frame of the demo script, and did some modifications, like 1)defining a different velocity field, 2)changing boundary condition and source term as a manufactured problem, 3)reducing to a stationary problem. I used the same SUPG stabilization term as in the demo, and expected a good solution.
However, the calculated result is not consistent with the exact solution. I can't figure out which part goes wrong in my script. It will be great if you give me a clue.
A short working script is attached here.

'''
Problem:
    b*grad(u) - div(c*grad(u)) = f  in domain
                             u = g  on boundary
    where b = [1, 1], c = 1, f = cos(x[0]) + cos(x[1]) + sin(x[0]) + sin(x[1]), and g = 'sin(x[0]) + sin(x[1])'
Exact solution;
    u = sin(x[0]) + sin(x[1])
'''
from dolfin import *

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

# Define velocity as vector
velocity = as_vector([1.0, 1.0])

c = 1.0

# Test and trial functions
u, v = TrialFunction(Q), TestFunction(Q)

# Manufactured right hand side
f = Expression('cos(x[0]) + cos(x[1]) + sin(x[0]) + sin(x[1])')

# Residual
r = dot(velocity, grad(u)) - c*div(grad(u)) - f

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

# Add SUPG stabilisation terms
vnorm = sqrt(dot(velocity, velocity))
h = CellSize(mesh)
delta = h/(2.0*vnorm)
F += delta*dot(velocity, grad(v))*r*dx

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

# Define boundary condition
def u0_boundary(x, on_boundary):
    return on_boundary
U0_expression = Expression('sin(x[0]) + sin(x[1])')
bc = DirichletBC(Q, U0_expression, u0_boundary)

# Solve
u = Function(Q)
solve(a == L, u, bc)

# Print out calcuated and exact solution
u_exact = interpolate(U0_expression, Q)
print 'x--', 'u------------ ', 'u_exact-------'
for i in xrange(0, 11):
    x = [0.5, i/10.0]
    print i/10.0, u(x), u_exact(x)

plot(u, interactive = True)

The calculated result and exact solution:
x--- u-------------------- u_exact-----------
0.0 0.479425538604 0.479425538604
0.1 0.513434804217 0.579258955251
0.2 0.561320051723 0.678094869399
0.3 0.620772943328 0.774945745266
0.4 0.689892906796 0.868843880913
0.5 0.767492771946 0.958851077208
0.6 0.853337272127 1.044068012
0.7 0.948370475345 1.12364322584
0.8 1.05496498288 1.1967816295
0.9 1.17717296021 1.26275244823
1.0 1.32089652341 1.32089652341

asked Apr 11, 2015 by Chao Zhang FEniCS User (1,180 points)
edited Apr 11, 2015 by Chao Zhang

1 Answer

+2 votes
 
Best answer

Hi, the form without stabilization is missing the forcing term. You should have

F = v*dot(velocity, grad(u))*dx + c*dot(grad(v), grad(u))*dx - f*v*dx

The absence of f term in the form from demo has no effect since f = Constant(0.) there.

answered Apr 11, 2015 by MiroK FEniCS Expert (80,920 points)
selected Apr 11, 2015 by Chao Zhang

Thank you very much for your time. I'm really sorry I missed that part.

...