Search code examples
pythonnumpyscipymayavifipy

Using FiPy and Mayavi to solve the diffusion equation in 3D

I'm interested in solving,

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

The following is working, but I have a few questions:

  1. Is it possible to increase performance with FiPy? I feel like the nx, ny, nz bins are very small here, despite a long computation time. I don't understand why the arrays X, Y, and Z are so large.
  2. Notice in the first frame, we are zoomed in. How can I force the extents to automatically be [0..nx, 0..ny, 0..nz] in all plots?
  3. Data for the first frame is a sphere of points with values 1.0 surrounded by 0.0. Why does there appear to be a gradient? Is Mayavi interpolating? If so, how can I disable this?

Code:

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()

Time elapsed: 68.2 seconds

Zeroth frame

First frame

Second frame

Third frame

Solution

    1. It's hard to tell from your question, but in the course of diagnosing things, I discovered that the LinearLUSolver scales very poorly as the dimension of the problem increases (see https://github.com/usnistgov/fipy/issues/474).

      For this symmetric problem, PySparse should use the PCG solver and Trilinos should use GMRES. If you didn't install either of these, then you'll get the SciPy sparse solvers, which defaults to LU (I don't know why; something for us to look into), and things will be really slow in 3D. Try adding solver=LinearGMRESSolver() to your eq.solve(...) statement.

      As far as the size of X, Y, and Z, you've declared a 30*30*30 cube of cells, so each of the cell center coordinate vectors will be 27000 elements long. Did you have a different expectation for cellCenters?

    2. I suggest you subclass our MayaviDaemon class, or at least look at how it sets up the display in Mayavi. In short, we set a data_set_clipper to the desired bounds.

    3. I don't know.