How To Couple Advection Diffusion Reaction Pdes With Fipy

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

Solution 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