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.

  • $\begingroup$ 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. $\endgroup$
    – Kbzon
    Commented Sep 22, 2015 at 15:01
  • $\begingroup$ @SantiagoLópezCastaño -- err, the Burgers equation is nonlinear, but the one in question here is still an advection diffusion equation. $\endgroup$ Commented Sep 23, 2015 at 14:41
  • $\begingroup$ That said, this is a question for the fipy project's mailing lists and forums. $\endgroup$ Commented Sep 23, 2015 at 14:42
  • 1
    $\begingroup$ @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? $\endgroup$ Commented Sep 23, 2015 at 15:36
  • $\begingroup$ @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. $\endgroup$
    – Paul
    Commented Sep 24, 2015 at 2:13

1 Answer 1


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) + \
# 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?")

  • $\begingroup$ 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. $\endgroup$ Commented Oct 9, 2015 at 15:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.