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
