No flux boundary condition in 2D-diffusion problem

Asked by Kyle Wedgwood

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


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,


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)


End code

Question information

English Edit question
FEniCS Project Edit question
No assignee Edit question
Last query:
Last reply:
Revision history for this message
N.A. Borggren (nborggren) said :

       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.

Revision history for this message
Kyle Wedgwood (pmxkw2) said :

Hi Nathan,

Thanks for your reply.

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.


Revision history for this message
N.A. Borggren (nborggren) said :

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.


Can you help with this problem?

Provide an answer of your own, or ask Kyle Wedgwood for more information if necessary.

To post a message you must log in.