python - 在 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 方法，您正在为 transient 解提供初始条件(或直接稳态的初始猜测)状态解决)而不是实际来源。所以正确的方法是你的第二种方法，你实际上直接将源/汇项添加到等式中:

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

此外，要尝试直接稳态求解，您可以在不使用 TransientTerm 的情况下编写方程式，

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

然后使用 eq.solve() 求解，这将求解码合的 DiffusionTerm 和源/汇。

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

dc/dt = R

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

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

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

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