Code related to FiPy mailing list question: http://article.gmane.org/gmane.comp.python.fipy/3551

convection.py

import matplotlib.pyplot as plt
from matplotlib import colors
from numpy.random import normal
from matplotlib.mlab import bivariate_normal
from fipy import *
from fipy import numerix
import numpy as np

fig = plt.figure
ax1 = plt.subplot((111))

# mesh config
l = 8.
nx = 201.
ny = nx
dx = l/nx
dy = dx
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy) + [[-l/2]]

convection = VanLeerConvectionTerm

# initial condition
x, y = mesh.getCellCenters()
z = bivariate_normal(x, y, .5, .5, 0., 0.)

# variable config
phi = CellVariable(mesh=mesh, value=z)

# boundary conditions
facesTopRight = ((mesh.getFacesRight()) 
                | (mesh.getFacesTop()))
facesBottomLeft = ((mesh.getFacesLeft())
                | (mesh.getFacesBottom()))

BCs = (FixedFlux(faces=facesTopRight, value=0.),
       FixedFlux(faces=facesBottomLeft, value=0.))

# viewer config
viewer1 = MatplotlibViewer(vars=phi,
                           datamin = 0.0, 
                           limits={'xmin': -l/2, 'xmax': l/2, 'ymin': -l/2, 'ymax': l/2}, 
                           cmap = plt.cm.Spectral,
                           axes=ax1)

# parameters                
D = 1.0
c = 1.0
b = ((c,), (c,))

alpha = 2./(dx * 0.5)

pos = FaceVariable(mesh=mesh, value=mesh.getFaceCenters(), rank=1)

xv, yv = mesh.getFaceCenters()
r = numerix.sqrt((xv-1.)**2+(yv-1.)**2)

# define equation 
faceVelocity = b*numerix.tanh(alpha*r)*pos/r
v_max = np.max(np.array(abs(faceVelocity)))
eq3 = TransientTerm() == DiffusionTerm(coeff=D) + convection(coeff=faceVelocity)

# solver config
CFL = 0.1
dt = CFL * dx / v_max
steps = 1000

for step in range(steps):
    print 'step',step
    eq3.solve(var=phi,
              boundaryConditions = BCs,
              dt=dt)
    if step % 10 == 0:
        viewer1.plot()