How to Couple Advection Diffusion Reaction PDEs with FiPy

Question

I tried to solve 1D coupled PDEs for an advection-diffusion-reaction problem with the Matlab function Pdepe (https://www.mathworks.com/help/matlab/ref/pdepe.html). This function is not working properly in my case of a high advection term as compared to the diffusion term. Therefore, I searched and found this option of using the Python library FiPy to solve my PDEs system. My initial conditions are u1=1 for 4*L/10

My coupled equations are of the following form:

du1/dt = d/dx(D1 * du1/dx) + g * x * du1/dx - mu1 * u1 / (K + u1) * u2

du2/dt = d/dx(D2 * du2/dx) + g * x * du2/dx + mu2 * u1 / (K + u1) * u2

I tried to write it by combining the FiPy examples (examples.convection.exponential1DSource.mesh1D, examples.levelSet.advection.mesh1D, examples.cahnHilliard.mesh2DCoupled).

The following lines are not a working example but my first attempt to write the code. This is my first use of FiPy (out of the tests and examples of the documentation), so this might seem to miss the point completely for the regular users.

from fipy import *

g = 0.66
L = 10.
nx = 1000
mu1 = 1.
mu2 = 1.
K = 1.
D1 = 1.
D2 = 1.

mesh = Grid1D(dx=L / 1000, nx=nx)

x = mesh.cellCenters[0]
convCoeff = g*(x-L/2)

u10 = 4*L/10 < x < 6*L/10
u20 = 1.

u1 = CellVariable(name="u1", mesh=mesh, value=u10)
u2 = CellVariable(name="u2", mesh=mesh, value=u20)

## Neumann boundary conditions
u1.faceGrad.constrain(0., where=mesh.facesLeft)
u1.faceGrad.constrain(0., where=mesh.facesRight)
u2.faceGrad.constrain(0., where=mesh.facesLeft)
u2.faceGrad.constrain(0., where=mesh.facesRight)

sourceCoeff1 = -1*mu1*u1/(K+u1)*u2
sourceCoeff2 = 1*mu2*u1/(K+u1)*u2

eq11 = (TransientTerm(var=u1) == DiffusionTerm(coeff=D1, var=u1) + ConvectionTerm(coeff=convCoeff))
eq21 = (TransientTerm(var=u2) == DiffusionTerm(coeff=D2, var=u2) + ConvectionTerm(coeff=convCoeff))

eq12 = ImplicitSourceTerm(coeff=sourceCoeff1, var=u1)
eq22 = ImplicitSourceTerm(coeff=sourceCoeff2, var=u2)

eq1 = eq11 & eq12
eq2 = eq21 & eq22

eqn = eq1 & eq2
vi = Viewer((u1, u2))

for t in range(100):
    u1.updateOld()
    u2.updateOld()
    eqn.solve(dt=1.e-3)
    vi.plot()

Thank you for any suggestion or correction. If you happen to know a good tutorial for this specific kind of problem, I would be happy to read it, since I did not find anything better than the examples in the FiPy documentation.

Answer 1

Several issues:

• python chained comparisons do not work in numpy and therefore do not work in FiPy. So, write

  u10 = (4*L/10 < x) & (x < 6*L/10)
  

  Further, this makes u10 a field of Booleans, which confuses FiPy, so write

  u10 = ((4*L/10 < x) & (x < 6*L/10)) * 1.
  

  or, better yet, write

  u1 = CellVariable(name="u1", mesh=mesh, value=0., hasOld=True)
  u2 = CellVariable(name="u2", mesh=mesh, value=1., hasOld=True)
  u1.setValue(1., where=(4*L/10 < x) & (x < 6*L/10))
  

• ConvectionTerm takes a vector coefficient. One way to get this is

  convCoeff = g*(x-L/2) * [[1.]]
  

  which represents a 1D rank-1 variable
• If you declare which Variable a Term applies to, you must do it for all Terms, so write, e.g.,

  ConvectionTerm(coeff=convCoeff, var=u1)
  

• ConvectionTerm(coeff=g*x, var=u1) does not represent g * x * du1/dx. It represents d(g * x * u1)/dx. So, I believe you'll want

  ConvectionTerm(coeff=convCoeff, var=u1) - ImplicitSourceTerm(coeff=g, var=u1)
  

• ImplicitSourceTerm(coeff=sourceCoeff1, var=u1 does not represent -1*mu1*u1/(K+u1)*u2, rather it represents -1*mu1*u1/(K+u1)*u2*u1. So, for best coupling between equations, write

  sourceCoeff1 = -mu1*u1/(K+u1)
  sourceCoeff2 = mu2*u2/(K+u1)
  
  ... ImplicitSourceTerm(coeff=sourceCoeff1, var=u2) ...
  ... ImplicitSourceTerm(coeff=sourceCoeff2, var=u1) ...
  

• As pointed out by @wd15 in the comments, you are declaring four equations for two unknowns. & does not mean "add two equations together" (which can be accomplished with +), rather it means "solve these two equations simultaneously". So, write

  sourceCoeff1 = mu1*u1/(K+u1)
  sourceCoeff2 = mu2*u2/(K+u1)
  
  eq1 = (TransientTerm(var=u1) 
         == DiffusionTerm(coeff=D1, var=u1) 
         + ConvectionTerm(coeff=convCoeff, var=u1) 
         - ImplicitSourceTerm(coeff=g, var=u1) 
         - ImplicitSourceTerm(coeff=sourceCoeff1, var=u2))
  eq2 = (TransientTerm(var=u2) 
         == DiffusionTerm(coeff=D2, var=u2) 
         + ConvectionTerm(coeff=convCoeff, var=u2) 
         - ImplicitSourceTerm(coeff=g, var=u2) 
         + ImplicitSourceTerm(coeff=sourceCoeff2, var=u1))
  
  eqn = eq1 & eq2
  

• A CellVariable must be declared with hasOld=True in order to call updateOld(), so

  u1 = CellVariable(name="u1", mesh=mesh, value=u10, hasOld=True)
  u2 = CellVariable(name="u2", mesh=mesh, value=u20, hasOld=True)
  

Full code that seems to work is here