Diffusion with space dependent drift in Fipy

I need to solve a diffusion equation in periodic boundary conditions using fipy but I would like to have a drift term that depends on the position so like this: $$\partial_t u(x,t) = \partial_x^2 u(x,t) + \partial_x(\omega(x)u(x,t)).$$

I am new in fipy and I am able to solve diffusion with diffusion coefficient spatially dependent but I am having troubles with this problem.

Does someone know how to manage? Or maybe if you know an easier way to solve numerically this equation.

• Your equation is no longer a diffusion equation, but a Burgers'-type law (of parabolic nature). Depending on the relative importance of the convective term, you may have regions with very steep gradients which are better solved using TVD or Godunov Schemes. – Kbzon Sep 22 '15 at 15:01
• @SantiagoLópezCastaño -- err, the Burgers equation is nonlinear, but the one in question here is still an advection diffusion equation. – Wolfgang Bangerth Sep 23 '15 at 14:41
• That said, this is a question for the fipy project's mailing lists and forums. – Wolfgang Bangerth Sep 23 '15 at 14:42
• @Paul, I think we are a bit too trigger-happy with closing software questions -- this one would have been on topic, I'd say ("how do I solve this scientific problem with this software", not "I get this error message, what does it mean"), even if it could have been improved. Maybe it's time for another discussion on meta about where we want to draw the line? – Christian Clason Sep 23 '15 at 15:36
• @ChristianClason: Upon re-examination, I agree with you. The question needs to be improved a bit, but should otherwise be on-topic. Indeed, the topicality of software-specific questions remains contentious and certainly merits further discussion. – Paul Sep 24 '15 at 2:13

A variable convection term works in the same way as a variable diffusion term:1 for (say) $\omega(x) = 1+x^2$ (and Dirichlet conditions for simplicity), you would write

from fipy import *
# generate mesh
nx = 50
dx = 1.
mesh = Grid1D(nx=nx, dx=dx)
x = mesh.cellCenters[0]
# solution phi on mesh with boundary conditions phi(0)=0, phi(1)=1
phi = CellVariable(name="solution variable", mesh=mesh, value=0.)
phi.constrain(0., mesh.facesRight)
phi.constrain(1., mesh.facesLeft)
# define pde phi_t = phi_xx + ((1+x^2)*phi)_x
D = 1.
omega = 1.+x**2
eqX = TransientTerm() == ExplicitDiffusionTerm(coeff=D) + \
ConvectionTerm(CellVariable(mesh=mesh,value=[omega]))
# time steps
timeStepDuration = 0.9 * dx**2 / (2 * D)
steps = 100
t = timeStepDuration * steps
# solve
for step in range(steps):
eqX.solve(var=phi, dt=timeStepDuration)


1. see http://ctcms.nist.gov/fipy/documentation/FAQ.html (under "What if my term involves the dependent variable, but not where FiPy puts it?")

• Thanks! I guess it works, but I am not able to plot the solution. I would like to plot the solution u(x,t) vs x for different time points. – user3585292 Oct 9 '15 at 15:54