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)
