FiPy 中使用 ConvectionTerm() 或 DiffusionTerm() 的数学运算

Question

我目前正在学习如何使用 FiPy，并最终想用它来解决一些生物学问题。我一直在尝试实现以下描述真菌菌丝生长的 PDE 系统：

PDE system

我可以毫无问题地实现这个系统，只要我忽略了si随时间的变化取决于第四个方程dsi中在空间abs(dp/dx)上的p的绝对变化 /dt。这是我没有 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()

现在，我试着简单地写

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))))

其中（预期）给出错误：“abs() 的操作数类型错误：'PowerLawConvectionTerm'”

我尝试通过添加另一个 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)

然后给出错误“ValueError：具有多个元素的数组的真值不明确。使用 a.any() 或 a.all()”，这可能是由于 eqdpdx 不是导函数？

我的问题是：是否可以使用 ConvectionTerm 执行数学运算？我能以某种方式表达 p 的绝对变化吗？ 而且，我如何表达像 eqdpdx（或任何动态参数）这样随空间变化的非导数函数？ 我查看了 FiPy 手册，但找不到解决方案 - 如果我的问题微不足道或已在其他地方得到解答，我很抱歉。

Answer 1

了解 FiPy Terms 是什么很重要。它们是 PDE 部分的人类可读表达式，可以离散化为线性代数。如果您对该线性代数应用一些潜在的非线性函数，那么它就不再是线性代数了。不支持这种用法，我什至不确定它是如何实现的。

幸运的是，您不需要它。 dp/dx 不是ConvectionTerm，它是一个渐变。我会把这个词写成

- ImplicitSourceTerm(coeff=c4*Da*m*p.grad.mag, var=si)

顺便说一句，一维方程是魔鬼的杰作。他们每次都会让你误入歧途。您可以像我们在手册中那样使用 nabla 表示法，也可以使用 Einstein 表示法，但请始终牢记您的表达是标量还是矢量（或张量或...）。