使用FiPy和Mayavi求解3D扩散方程
Using FiPy and Mayavi to solve the diffusion equation in 3D
我有兴趣解决,
\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 秒
很难从你的问题中判断出来,但在诊断过程中,我发现 LinearLUSolver
缩放 非常差 作为问题的维度增加(参见 https://github.com/usnistgov/fipy/issues/474)。
对于这个对称问题,PySparse 应该使用 PCG 求解器,Trilinos 应该使用 GMRES。如果你没有安装其中任何一个,那么你将获得 SciPy 稀疏求解器,默认为 LU(我不知道为什么;我们需要研究的东西),事情会很慢在 3D 中。尝试将 solver=LinearGMRESSolver()
添加到您的 eq.solve(...)
语句中。
就 X、Y 和 Z 的大小而言,您已经声明了一个 30*30*30 的单元格立方体,因此每个单元格中心坐标向量的长度为 27000 个元素。您对 cellCenters
有不同的期望吗?
我建议你subclass我们的MayaviDaemonclass,或者至少看看它是如何在Mayavi中设置显示的。简而言之,我们将 data_set_clipper
设置为所需的范围。
我不知道。
我有兴趣解决,
\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 秒
很难从你的问题中判断出来,但在诊断过程中,我发现
LinearLUSolver
缩放 非常差 作为问题的维度增加(参见 https://github.com/usnistgov/fipy/issues/474)。对于这个对称问题,PySparse 应该使用 PCG 求解器,Trilinos 应该使用 GMRES。如果你没有安装其中任何一个,那么你将获得 SciPy 稀疏求解器,默认为 LU(我不知道为什么;我们需要研究的东西),事情会很慢在 3D 中。尝试将
solver=LinearGMRESSolver()
添加到您的eq.solve(...)
语句中。就 X、Y 和 Z 的大小而言,您已经声明了一个 30*30*30 的单元格立方体,因此每个单元格中心坐标向量的长度为 27000 个元素。您对
cellCenters
有不同的期望吗?我建议你subclass我们的MayaviDaemonclass,或者至少看看它是如何在Mayavi中设置显示的。简而言之,我们将
data_set_clipper
设置为所需的范围。我不知道。