使用 FiPy 和 Mayavi 求解 3D 中的扩散方程答案

【问题标题】：Using FiPy and Mayavi to solve the diffusion equation in 3D使用 FiPy 和 Mayavi 求解 3D 中的扩散方程
【发布时间】：2015-11-17 18:54:40
【问题描述】：

我有兴趣解决，

\frac{\delta \phi}{\delta t} - D \nabla^2 \phi - \alpha \phi - \gamma \phi = 0

以下是有效的，但我有几个问题：

是否可以使用 FiPy 提高性能？我觉得 nx, ny, nz 这里的垃圾箱非常小，尽管计算时间很长。 ~~我不明白为什么数组X, Y, and Z 这么大。~~
注意在第一帧中，我们被放大了。如何强制所有地块中的范围自动为[0..nx, 0..ny, 0..nz]？
第一帧的数据是一个球体，其值为1.0，被0.0 包围。为什么会出现渐变？ Mayavi 是在插值吗？如果是这样，我该如何禁用它？

代码：

from fipy import *
import mayavi.mlab as mlab
import numpy as np
import time

# Spatial parameters
nx = ny = nz = 30  # bins
dx = dy = dz = 1  # Must this be an integer?
L = nx * dx

# Diffusion and time step
D = 1.
dt = 10.0 * dx**2 / (2. * D)
steps = 4

# Initial value and radius of concentration
phi0 = 1.0
r = 3.0

# Rates
alpha = 1.0  # Source coeficcient
gamma = .01  # Sink coeficcient

mesh = Grid3D(nx=nx, ny=ny, nz=nz, dx=dx, dy=dy, dz=dz)
X, Y, Z = mesh.cellCenters  # These are large arrays
phi = CellVariable(mesh=mesh, name=r"$\phi$", value=0.)

src = phi * alpha  # Source term (zeroth order reaction)
degr = -gamma * phi  # Sink term (degredation)

eq = TransientTerm() == DiffusionTerm(D) + src + degr

# Initial concentration is a sphere located in the center of a bounded cube
phi.setValue(1.0, where=( ((X-nx/2))**2 + (Y-ny/2)**2 + (Z-nz/2)**2 < r**2) )

# Solve
start_time = time.time()
results = [phi.getNumericValue().copy()]
for step in range(steps):
    eq.solve(var=phi, dt=dt)
    results.append(phi.getNumericValue().copy())
print 'Time elapsed:', time.time() - start_time

# Plot
for i, res in enumerate(results):
    fig = mlab.figure()

    res = res.reshape(nx, ny, nz)
    mlab.contour3d(res, opacity=.3, vmin=0, vmax=1, contours=100, transparent=True, extent=[0, 10, 0, 10, 0, 10])

    mlab.colorbar()
    mlab.savefig('diffusion3d_%i.png'%(i+1))
    mlab.close()

经过时间：68.2 秒

【问题讨论】：

标签： python numpy scipy mayavi fipy

【解决方案1】：

很难从您的问题中看出，但在诊断过程中，我发现LinearLUSolver 随着问题规模的增加非常差（请参阅https://github.com/usnistgov/fipy/issues/474）。

对于这个对称问题，PySparse 应该使用 PCG 求解器，而 Trilinos 应该使用 GMRES。如果您没有安装其中任何一个，那么您将获得默认为 LU 的 SciPy 稀疏求解器（我不知道为什么；我们需要研究一下），并且在 3D 中事情会非常慢。尝试将solver=LinearGMRESSolver() 添加到您的eq.solve(...) 语句中。

就 X、Y 和 Z 的大小而言，您已经声明了一个 30*30*30 的单元格立方体，因此每个单元格中心坐标向量的长度为 27000 个元素。你对cellCenters有什么不同的期待吗？
我建议您将我们的MayaviDaemon 类子类化，或者至少看看它是如何在 Mayavi 中设置显示的。简而言之，我们将data_set_clipper 设置为所需的范围。
我不知道。

【讨论】：