2

我正在使用 FiPy 解决受生物学启发的问题。

本质上，我想代表一个 2D 平面，在不同的点我有源和汇。源以固定的速率发射底物（不同的源可以有不同的速率），而汇以固定的速率消耗底物（不同的汇可以有不同的速率）。我的代码：

import numpy.matlib
from fipy import CellVariable, Grid2D, Viewer, TransientTerm, DiffusionTerm, ImplicitSourceTerm, ExplicitDiffusionTerm
from fipy.tools import numerix
from time import *

nx = 10
ny = nx
dx = 1.
dy = dx
L = dx*nx
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

arr_grid = numerix.array((
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,),'d')

arr_source = numerix.array((
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0.5,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,),'d')

arr_sink = numerix.array((
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0.5,),'d')

source = CellVariable(mesh=mesh, value = arr_source)
sink = CellVariable(mesh=mesh, value = arr_sink)
phi = CellVariable(name = "solution variable", mesh = mesh, value = arr_grid)
X,Y = mesh.cellCenters
phi.setValue(3.0, where=(X < 2.0) & (X > 1.0))
phi.setValue(-1.0, where=(X < 6.0) & (X > 5.0))
D = 1.
eq = TransientTerm() == DiffusionTerm(coeff=D)



viewer = Viewer(vars=phi, datamin=0., datamax=1.)

steadyState = False

if(steadyState):
    print("SteadyState")
    DiffusionTerm().solve(var=phi)
    viewer.plot()
    sleep(20)
else:
    print("ByStep")
    timeStepDuration = 10 * 0.9 * dx**2 / (2 * D)
    steps = 500
    for step in range(steps):
        print(step)
        eq.solve(var=phi,
        dt=timeStepDuration)
        if __name__ == '__main__':
            viewer.plot()

这很好用，但 FiPy 将来源视为“不可再生”，最终我在整个空间中得到了预期的均匀浓度。另一种方法是删除：

X,Y = mesh.cellCenters
phi.setValue(3.0, where=(X < 2.0) & (X > 1.0))
phi.setValue(-1.0, where=(X < 6.0) & (X > 5.0))

并将等式更改为：

eq = TransientTerm() == DiffusionTerm(coeff=D) + source - sink

鉴于源和汇永远不会改变，这提供了“无限”的源和汇。但是，当我尝试使用求解稳态时

eq = TransientTerm() == DiffusionTerm(coeff=D) + source - sink

我得到：

C:\Python27\python.exe C:/Users/dario_000/PycharmProjects/mesa-test/mesher.py
SteadyState
C:\Python27\lib\site-packages\fipy-3.1.dev134+g64f7866-py2.7.egg\fipy\solvers\scipy\linearLUSolver.py:71: RuntimeWarning: invalid value encountered in double_scalars
  if (numerix.sqrt(numerix.sum(errorVector**2)) / error0)  <= self.tolerance:

而且方程没有解。但是，如果我再次使用“逐步”解决它：

eq = TransientTerm() == DiffusionTerm(coeff=D) + source - sink

我得到了一张类似于我所期待的漂亮图片：

在此处输入图像描述

关于如何在不同的空间位置指定具有不同排放/消耗率的源/汇的初始设置的任何建议，以便我可以获得稳态解决方案？

谢谢！

4

1 回答 1

4

请注意，正如 wd15 在评论中提到的，邮件列表上有更深入的讨论和跟进。

首先，初始条件和来源都可以使用where逻辑来制作。

source = CellVariable(mesh=mesh, value = arr_source, where=(2 < X) & (X < 3))

这避免了数组的显式构造。

其次，初始条件和来源之间存在关键区别。在原始公式中，您SetValue在字段变量 上使用方法phi，您为瞬态解（或直接稳态解的初始猜测）提供初始条件，而不是实际源。所以正确的方法是你的第二种方法，你实际上将源/汇项直接添加到方程中：

eq = TransientTerm() == DiffusionTerm(coeff=D) + source - sink

此外，要尝试直接稳定求解，您可以编写不带TransientTerm,

eq = 0 == DiffusionTerm(coeff=D) + source - sink

然后解决 using eq.solve()，这将解决组合DiffusionTerm和源/汇。

然而，直接稳定的方法值得几句谨慎的话。首先，没有一种很好的数值方法可以直接计算一般系统的稳态解。通常，您最好的选择实际上是通过从某个初始条件到达到稳态的时间步长来求解瞬态方程，因为这实际上可能是求解稳态曲线的最稳健的算法。在某些情况下，你甚至可以用一个大的时间步长来做到这一点，就像在“完全隐式解决方案并非没有陷阱”部分开始的那样. 其次，您确定您的系统承认稳定状态吗？您有无通量边界条件（暗示因为您没有指定任何其他边界条件），但您有内部源/汇。除非这些源/汇以完全相同的速率产生/消耗场变量，否则您将在系统中获得净积累。一个带有R= 常量、统一和非零的简单示例：

dc/dt = R

没有通量边界条件是一个不允许任何稳态的方程。添加扩散项不会改变这一点。如果您要添加任何狄利克雷（指定值）边界条件，这将实现稳定的解决方案，因为系统内的净生产/消耗可以通过系统边界离开/进入。可以在此处找到有关边界条件以及如何应用它们的讨论。

最后，还有一点需要注意。所写的源/汇项是“零阶”，如果汇项足够大和/或离源足够远，这将导致例如浓度为负。如果发生这种情况，模型显然需要修改，例如，将接收器设为一阶，

eq = TransientTerm() == DiffusionTerm(coeff=D) + source - sink*phi

这将确保当 phi 变为零时汇“关闭”，但这也可以通过修改以耦合到其他场变量，如局部细胞密度。

于 2016-04-05T20:20:07.590 回答