使用 Python 中的源解决 PDE

Question

我正在使用 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

我得到了一张与我预期相似的漂亮图片：

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

谢谢！

Answer 1

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

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

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

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

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

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

另外，要尝试直接稳定求解，您可以编写不带TransientTerm 的方程，

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

然后使用eq.solve() 求解，这将求解组合的DiffusionTerm 和源/接收器。

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

dc/dt = R

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

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

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

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