Hi all,
I am quite new to FiPy, which I have been using to simulate a reaction-
diffusion model of bacterial colony growth:
# Simulation parameters
d = 0.12
a0 = 1.; a1 = 1./2400; a2 = 1./120;
v0 = 0.071
# Source terms
Suu = uu*vv - uu/((1+uu/a1)*(1+vv/a2))
Svv = -uu*vv
Sww = uu/((1+uu/a1)*(1+vv/a2))
# Equations
eqXuu = TransientTerm() - ImplicitDiffusionTerm(coeff=d) - Suu
eqXvv = TransientTerm() - ImplicitDiffusionTerm(coeff=1.) - Svv
eqXww = TransientTerm() - Sww
There are three PDEs, two of which have linear diffusion terms, and
the equations are coupled via nonlinear reaction terms. I hope to be
able to use Fipy on a fairly large mesh to explore the pattern
formation possibilities of the system as certain parameters are
varied, but the simulations run significantly slower than I'd
expected, about 4-5 minutes per run (with --inline and the pysparse
solver on a 2008 Mac Pro) on a 150x150 square mesh; 600x600 is more
like what I'd like to be able to use.
Before using FiPy, I'd implemented the system in Matlab, just using
finite differences with a 9 point stencil. I haven't gotten much, if
any, speed up using FiPy, but I think I'm not using FiPy's full
capabilities, or I've not implemented the model in FiPy as efficiently
as I could have. The script I'm using is below. I suspect that the
looping I'm doing for the sweeps of each variable in succession is not
the right way to go about things. I've played around with sweeps and
time step sizes, but haven't gotten much improvement.
If the looping isn't the problem, should I try moving to trilinos? Or
should I not expect particularly fast solution of this kind of system?
Any suggestions for better performance would be welcome.
Thanks,
Erik
from fipy import *
from time import clock
# Setting up the mesh
nx = 150
ny = nx
dx = 1.0
dy = dx
L = dx*nx
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
# Simulation parameters
d = 0.12
a0 = 1.; a1 = 1./2400; a2 = 1./120;
v0 = 0.071
# Simulation variables
uu = CellVariable(name="active bacteria",
mesh=mesh,
value=0.0, hasOld=1)
ww = CellVariable(name="inactive bacteria",
mesh=mesh,
value=0.0, hasOld=1)
vv = CellVariable(name="nutrient concentration",
mesh=mesh,
value=v0, hasOld=1)
# Source terms
Suu = uu*vv - uu/((1+uu/a1)*(1+vv/a2))
Svv = -uu*vv
Sww = uu/((1+uu/a1)*(1+vv/a2))
# Equations
eqXuu = TransientTerm() - ImplicitDiffusionTerm(coeff=d) - Suu
#ImplicitSourceTerm(Suu)
eqXvv = TransientTerm() - ImplicitDiffusionTerm(coeff=1.) - Svv
#ImplicitSourceTerm(Svv)
eqXww = TransientTerm() - Sww #ImplicitSourceTerm(Sww)
# Initial conditions
x,y = mesh.getCellCenters()
radius = 10*dx
C = (nx*dx/2, ny*dy/2)
C =(0,0)
uu.setValue(UniformNoiseVariable(mesh=mesh, minimum=0.09,
maximum=0.115))
uu.setValue(0.0,where=(0.05*(x-C[0])**2 + (y-C[1])**2) > radius**2)
eqns = ((uu, eqXuu), (vv, eqXvv), (ww, eqXww))
if __name__ == '__main__':
# Do plotting if executing script from commandline
Uview = uu
Uview.setName('Active')
Uviewer = Viewer(vars=uu, datamin=0., datamax=0.1)
Vview = vv
Vview.setName('Nutrient')
Vviewer = Viewer(vars=Vview, datamin=0., datamax=0.1)
Wview = ww
Wview.setName('Inactive')
Wviewer = Viewer(vars=Wview, datamin=0., datamax=0.1)
dt = 2.0
steps = 1000
desiredResidual = 1e-6
# Loop through timesteps
tic = clock()
for step in range(steps):
print "Step %d"%(step)
# For each variable, equation pair, do sufficient sweeps to
# get residual to acceptable level. When all three variables
# have been taken care of, then we take the next time step
for var, eqn in eqns:
residual = 10
var.updateOld()
sweepcount = 0
while residual > desiredResidual and sweepcount < 300:
#print "Sweep %d"%(sweepcount+1)
residual = eqn.sweep(var=var, dt=dt)
sweepcount += 1
# Bail if residual remains too large
if residual > desiredResidual:
print "Unable to reduce residual!"
break
if residual > desiredResidual:
break
# Update the plots
if __name__ == '__main__':
Uviewer.plot()
Vviewer.plot()
Wviewer.plot()
#
toc = clock()
print "Finished in %.3f secs."%(toc-tic)
