2015-11-17 48 views
0

擴散方程我感興趣的解決,使用FiPy和Mayavi的解決3D

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

以下是工作,但我有幾個問題:

  1. 是否有可能通過FiPy提高性能?儘管計算時間很長,但我覺得這裏的nx, ny, nz箱子非常小。 我不明白爲什麼陣列 X, Y, and Z是如此之大。
  2. 請注意,在第一幀中,我們放大了。如何強制範圍在所有地塊中自動爲[0..nx, 0..ny, 0..nz]
  3. 第一個幀的數據是一個由0.0包圍的值爲1.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秒

Zeroth frame

First frame

Second frame

Third frame

回答

4
  1. 很難從你的問題告訴,但在診斷事物的過程中,我發現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有不同的期待嗎?

  2. 我建議你繼承我們的MayaviDaemon類,或者至少看看它如何在Mayavi中設置顯示。簡而言之,我們將data_set_clipper設置爲期望的範圍。

  3. 我不知道。