# 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

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

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)
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

## Question information

Language:
English Edit question
Status:
For:
FEniCS Project Edit question
Assignee:
No assignee Edit question
Last query:
 Revision history for this message N.A. Borggren (nborggren) said on 2013-03-05: #1

Hello,
To impose the no flux boundary condition, the term - 2*sigma*dt*v*(inner(dot(A,grad(u))),n)*ds needs to be set to zero (by removing the term). Have a look at the advection-diffusion demos, I think they would be very useful for your case.
Cheers,
Nathan

 Revision history for this message Kyle Wedgwood (pmxkw2) said on 2013-03-05: #2

Hi Nathan,

I've set the term in question to zero, and I my solution still blows up (even with a small step size)... I'm still very confused by why this keeps happening.

Kye

 Revision history for this message N.A. Borggren (nborggren) said on 2013-03-05: #3

Have a look at what the "midpoint solution" is trying to do in the advection-diffusion demo. Replacing your Q with (Q+Q1)/2 and collecting terms with Q in your "a" form and terms with "Q1" in the L form might work better then the forms you choose.
Also try an initial condition that steers clear of the boundary, a gaussian for example, in my experience constant initial conditions are difficult to implement. That and using the finest mesh and smallest time step your computer can support might be necessary. It is indeed a finicky equation, as evidenced by all the efforts for stabilization in the demo.

Nathan