I'm currently learning how to use FiPy and eventually want to use it to solve some biological problems. I have been trying to implement the following PDE system which describes hyphal growth of fungi:
I can implement the system without problems, as long as I ignore that the change of si over time depends on the absolute change of p over space abs(dp/dx) in the fourth equation dsi/dt. Here is my code without abs():
from fipy import *
##### produce mesh
nx= 100.
ny= nx
dx = 1/100
dy = dx
L = nx*dx
mesh = Grid2D(nx=nx,ny=ny,dx=dx,dy=dy)
x,y = mesh.cellCenters
##### parameters
b=1e+7 # branching rate (no. branches * cm^-1 * hyphae * day^-1 * (mol glucose)^-1)
f = 10. # fusion rate (no. fusions in cm * day^-1)
v = 1e+5 # cm * mol^-1 * day^-1
d = 0.5 # day^-1
r = 0. # degradation of hyphae
c1 = 9e+2 # cm * mol^-1 * day^-1
c2 = 1e-7 # mol * cm^-1
c3 = 1e+3 # cm * mol^-1 * day^-1, c1/c3 = 90% efficiency
c4 = 1e-8 # cm^-1
De = 1e-3 # diffusion of external substrate (cm² * s^-1)
Di = 1e-2 # diffusion of internal substrate (cm³ * day^-1)
Da = 0.
##### define state variables:
m = CellVariable(name="m",mesh=mesh,hasOld=True,value=0.)
mi = CellVariable(name="mi",mesh=mesh,hasOld=True,value=0.)
p = CellVariable(name="p",mesh=mesh,hasOld=True,value=0.)
si = CellVariable(name="si",mesh=mesh,hasOld=True,value=0.)
se = CellVariable(name="se",mesh=mesh,hasOld=True,value=0.)
##### differential equations
eqm = (TransientTerm(var=m)== si*v*p - ImplicitSourceTerm(var=m,coeff=d))
eqmi = (TransientTerm(var=mi)==m*d - ImplicitSourceTerm(var=mi,coeff=r))
eqp = (TransientTerm(var=p)== - ConvectionTerm(var=p,coeff=[[v]]*si) + b*si*m - ImplicitSourceTerm(var=p,coeff=f*m))
eqsi = (TransientTerm(var=si)== DiffusionTerm(var=si,coeff=Di*m)- DiffusionTerm(var=p,coeff=Da*m*si) + ImplicitSourceTerm(var=si,coeff=c1*m*se) - ImplicitSourceTerm(var=si,coeff=c2*v*p) - ConvectionTerm(var=p,coeff=[[c4*Da]]*(m*si)))
eqse = (TransientTerm(var=se) == DiffusionTerm(var=se,coeff=De) - ImplicitSourceTerm(var=se,coeff=c3*m*si))
# initial values
m0 = 100.
mi0 = 0.
p0 = 500.
si0 = 1e-5
se0 = 3e-5
r = 3. # radius of initial plug
m.setValue(m0,where=((x>(L/2-(r*dx)))&(x<(L/2+(r*dx)))&(y>(L/2-(r*dy)))&(y<(L/2+(r*dy)))))
mi.setValue(mi0,where=((x>(L/2-(r*dx)))&(x<(L/2+(r*dx)))&(y>(L/2-(r*dy)))&(y<(L/2+(r*dy)))))
p.setValue(p0,where=((x>(L/2-(r*dx)))&(x<(L/2+(r*dx)))&(y>(L/2-(r*dy)))&(y<(L/2+(r*dy)))))
si.setValue(si0,where=((x>(L/2-(r*dx)))&(x<(L/2+(r*dx)))&(y>(L/2-(r*dy)))&(y<(L/2+(r*dy)))))
se.setValue(se0)
# boundary conditions
#----------------
# simulation
eq = eqm & eqmi & eqp & eqsi & eqse
vi = Viewer((m))
from builtins import range
for t in range(100):
m.updateOld()
mi.updateOld()
p.updateOld()
si.updateOld()
se.updateOld()
eq.solve(dt=0.1)
print(t)
vi.plot()
Now, I tried to simply write
eqsi = (TransientTerm(var=si)== DiffusionTerm(var=si,coeff=Di*m)- DiffusionTerm(var=p,coeff=Da*m*si) + ImplicitSourceTerm(var=si,coeff=c1*m*se) - ImplicitSourceTerm(var=si,coeff=c2*v*p) - abs(ConvectionTerm(var=p,coeff=[[c4*Da]]*(m*si))))
which (expectedly) gives out an error: "bad operand type for abs(): 'PowerLawConvectionTerm'"
I tried to work my way around it by adding another CellVariable dpdx:
dpdx= CellVariable(name="dpdx",mesh=mesh,hasOld=True,value=0.)
eqdpdx= (dpdx == ConvectionTerm(var=p,coeff=[[1]]))
eqsi = (TransientTerm(var=si)== DiffusionTerm(var=si,coeff=Di*m)- DiffusionTerm(var=p,coeff=Da*m*si) + ImplicitSourceTerm(var=si,coeff=c1*m*se) - ImplicitSourceTerm(var=si,coeff=c2*v*p) - abs(dpdx)*c4*Da*m*si)
which then gives the error "ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()", which is probably caused by the fact that eqdpdx is not a derivative function?
My question is: Is it possible at all to perform mathematical operations with a ConvectionTerm? Can I somehow express the absolute change of p? And also, how can I express non derivative function that change over space like eqdpdx (or any dynamic parameters)? I looked through the FiPy manual and could not find a solution - I am sorry if my question is trivial or has been answered elsewhere.
It's important to understand what FiPy Terms are. They are human-readable expressions of parts of PDEs that can be discretized to linear algebra. If you apply some potentially non-linear function to that linear algebra, then it's not linear algebra anymore. This use is not supported and I'm not sure how it even could be.
Fortunately, you don't need it. dp/dx is not a ConvectionTerm, it's a gradient. I would write this term as
- ImplicitSourceTerm(coeff=c4*Da*m*p.grad.mag, var=si)
As a side note, 1D equations are the Devil's work. They will lead you astray every single time. You can use nabla notation like we do in the manual or you can use Einstein notation, but always keep clear in your head whether your expression is scalar or vector (or tensor or ...).