使用 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 将源视为 "non-renewable"，最终我在整个 space 中得到了预期的均匀集中。另一种方法是删除：

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

鉴于源和接收器从未更改，因此提供了 "infinite" 源和接收器。但是，当我尝试使用

求解稳态时

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:

而且方程没有解出来。但是，如果我再次使用 "by steps" 解决它：

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

我得到了一张与我期望的相似的漂亮照片：

关于如何在不同空间位置使用 sources/sinks 指定初始设置的任何建议，每个空间位置具有不同的 emission/consumption 速率，以便我可以获得稳态解决方案？

谢谢！

Answer 1

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

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

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

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

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

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

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

的方程式

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

然后使用 eq.solve() 求解，这将求解合并的 DiffusionTerm 和 sources/sinks.

但是，直接到稳定的方法值得提醒几句。首先，没有真正好的数值方法来直接计算一般系统的稳态解。通常，您最好的选择实际上是通过从某个初始条件到达到稳态的时间步长来求解瞬态方程，因为这实际上可能是求解稳态曲线的最稳健的算法。在某些情况下，您甚至可以用一个大的时间步来做到这一点，如开头部分 "fully implicit solutions are not without their pitfalls" here 中所示。其次，你确定你的系统承认稳定状态吗？你有无通量边界条件（暗示是因为你没有指定任何其他边界条件），但你有内部 sources/sinks。除非那些 sources/sinks produce/consume 字段变量以完全相同的速率变化，否则您将在系统中有净积累。 R = 常量、统一且非零的简单示例：

dc/dt = R

没有通量边界条件是一个不允许任何稳态的方程。添加扩散项不会改变这一点。如果您要添加任何 dirichlet（指定值）边界条件，这将实现稳定的解决方案，因为系统内的网络 production/consumption 可以通过系统边界 leave/enter。可以找到关于边界条件和如何应用它们的讨论 here。

最后，还有一点要注意。 source/sink 条款，如书面的那样是 "zero-order,"，这将导致，例如如果汇项足够大 and/or 离源足够远，浓度将变为负值。如果发生这种情况，模型显然需要修改，例如，使汇一阶，

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

这将确保汇 "turn off" 随着 phi 变为零，但这也可以通过修改以耦合到其他场变量，如局部细胞密度。

使用 Python 中的源求解 PDE

Solving PDE with sources in Python

python

fipy