Eric, Thanks for your interest in FiPy. I'll benchmark your code, but
before I do that, can you tell me which version of FiPy you are using
(trunk or version-2_1?) so we are comparing apples with apples. Thanks
On Wed, Sep 29, 2010 at 5:03 PM, Erik Sherwood <wes...@bu.edu> wrote:
>
> 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)
>
>
>
>
--
Daniel Wheeler
