在 FiPy 中使用 ConvectionTerm() 或 DiffusionTerm() 进行数学运算
Mathematical operations with ConvectionTerm() or DiffusionTerm() in FiPy
我目前正在学习如何使用FiPy，最终想用它来解决一些生物学问题。我一直在尝试实施以下描述真菌菌丝生长的 PDE 系统：
PDE system
我可以毫无问题地实施该系统，只要我忽略 si 随时间的变化取决于 p 相对于 space abs 的 绝对变化（ dp/dx)在第四个方程dsi/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（或任何动态参数）那样随 space 变化的非导数函数？
我查看了 FiPy 手册但找不到解决方案 - 如果我的问题微不足道或已在其他地方得到解答，我深表歉意。
了解 FiPy Terms 是什么很重要。它们是可离散化为线性代数的 PDE 部分的人类可读表达式。如果您将一些潜在的非线性函数应用于该线性代数，那么它就不再是线性代数了。不支持这种用途，我什至不确定它会怎样。
幸运的是，你不需要它。 dp/dx 不是 ConvectionTerm，而是渐变。我会把这个词写成
- ImplicitSourceTerm(coeff=c4*Da*m*p.grad.mag, var=si)
附带说明一下，一维方程是魔鬼的杰作。他们每次都会让你误入歧途。您可以像我们在手册中那样使用 nabla 表示法，也可以使用 Einstein 表示法，但请始终牢记您的表达式是标量还是向量（或张量或...）。
我目前正在学习如何使用FiPy，最终想用它来解决一些生物学问题。我一直在尝试实施以下描述真菌菌丝生长的 PDE 系统：
PDE system
我可以毫无问题地实施该系统，只要我忽略 si 随时间的变化取决于 p 相对于 space abs 的 绝对变化（ dp/dx)在第四个方程dsi/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（或任何动态参数）那样随 space 变化的非导数函数？ 我查看了 FiPy 手册但找不到解决方案 - 如果我的问题微不足道或已在其他地方得到解答，我深表歉意。
了解 FiPy Terms 是什么很重要。它们是可离散化为线性代数的 PDE 部分的人类可读表达式。如果您将一些潜在的非线性函数应用于该线性代数，那么它就不再是线性代数了。不支持这种用途，我什至不确定它会怎样。
幸运的是，你不需要它。 dp/dx 不是 ConvectionTerm，而是渐变。我会把这个词写成
- ImplicitSourceTerm(coeff=c4*Da*m*p.grad.mag, var=si)
附带说明一下，一维方程是魔鬼的杰作。他们每次都会让你误入歧途。您可以像我们在手册中那样使用 nabla 表示法，也可以使用 Einstein 表示法，但请始终牢记您的表达式是标量还是向量（或张量或...）。