Использование FiPy и MayaVi для решения уравнения диффузии в 3D

Question

я заинтересован в решении,

\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. Данные для первого кадра - сфера точек со значениями 1.0, окруженная 0.0. Почему появляется градиент? Является ли Маяви интерполяцией? Если да, то как я могу отключить это?

Код:

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 секунд

Answer 1

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. Я не знаю.