No flux boundary condition in 2D-diffusion problem

Hi all,

I'm completely new to FE methods, so this may be a very simple problem, but I can't get my head around the issue.

I'm looking to solve the following time-dependent problem in two spatial variables (x and y)

\partial u/\partial t = -\partial{ f u }/\partial x -\partial{g u}/\partial y + 2*\sigma *\partial^2{u}/\partial{x^2},

on a square domain of size R, subject to no-flux boundary conditions and the initial condition

u(x,y) = 1/R^2 in \Omega.

I begin by transforming the equation into the form

\partial u/\partial t = -C(x,y)u - B(x,y),grad(u) + 2*sigma*div(A(x,y)*grad(u))

where

A = [ [ 1 , 0 ],[0 0 ] ] (matrix)
B = [ f , g ] ( vector)
C = \partial f/\partial x + \partial g/\partial y (scalar)

For the time stepping, I use a backward Euler scheme. Upon applying this, multiplying by a test function v and integrating by parts, I get

v*u*dx + dt*v*C*u*dx + dt*v*inner(B,grad(u))*dx - 2*sigma*dt*v*(inner(dot(A,grad(u))),n)*ds + 2*sigma*inner(grad(v),dot(A,grad(u)))*dx = v*u1*dx

where dt is my time step and, u1 is the solution from the previous time step and n is my outwardly pointing normal.

Since I have no flux boundary conditions, I expect the boundary term to disappear.

When, I try to solve this system, the solution simply diverges to infinity.

As a test case, I solve the system where

f = 0.5*lambda*x - (0.5*lambda*c+omega)*y - 0.5*(x^2+y^2)*(x-c*y)
g = (0.5*lambda*c+omega)*x + 0.5*lambda*y - 0.5*(x^2+y^2)*(c*x+y)

for which I expect the solution after a long time to look like a ring, centered on the unit circle.

Please can someone help my identify where I've gone wrong. I attach the code I use to solve the problem below.

Many thanks,

Kyle

Begin code:

from dolfin import *
from ufl import exp
import numpy

# define parameters
lam = 2.0
omega = 1.0
c = 1.0
sigma = 0.5

# create mesh
R = 5.0
mesh = RectangleMesh(-R,-R,R,R,100,100)

V = FunctionSpace(mesh,"Lagrange",1)

# define initial condition
Q0 = Constant(1/(pow(R,2)))
Q1 = interpolate(Q0,V)

# time stuff
T = 100.0
dt = 0.001

# Define variational problem
Q = TrialFunction(V)
v = TestFunction(V)

# functions
A = Constant(((1.0,0.0),(0.0,0.0)))
B = Expression(('0.5*lam*x[0]-(0.5*lam*c+omega)*x[1]-0.5*lam*(pow(x[0],2)+pow(x[1],2))*(x[0]-c*x[1])','(0.5*lam*c+omega)*x[0]+0.5*lam*x[1]-0.5*lam*(pow(x[0],2)+pow(x[1],2))*(c*x[0]+x[1])'),lam=lam,omega=omega,c=c)
C = Expression('lam*(1-2*pow(x[0],2)-2*pow(x[1],2))',lam=lam)
a = v*Q*dx+dt*v*C*Q*dx+dt*v*inner(B,nabla_grad(Q))*dx+2*sigma*dt*inner(dot(A,nabla_grad(Q)),nabla_grad(v))*dx
L = v*Q1*dx

A = assemble(a) # assemble only once, before the time stepping
b = None # necessary for memory saving assemble call
Q = Function(V) # the unknown at a new time level
t = dt

# solve problem
while t <= T:
    print 'time =', t
    b = assemble(L, tensor=b)
s solve(A, Q.vector(), b)
    t+=dt
    Q1.assign(Q)
    plot(Q)

interactive()

End code