I'm running a workshop right now, so won't have time to look at this for a day or two and it sounds like Daniel is on it.
One thing that occurs to me is that you should be able to use an ImplicitSourceTerm for Suu and Svv, which might reduce the number of sweeps. How many sweeps are you doing? How did you pick desiredResidual? Generally, you want to compare the residual to the initial residual, rather than some arbitrary number. On Sep 30, 2010, at 12:55 PM, Erik Sherwood wrote: > > Thanks for taking a look at the code. I'm using FiPy 2.1. > > Erik > > On Sep 30, 2010, at 11:21 AM, Daniel Wheeler wrote: > >> >> 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 <[email protected]> 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 >> >> >> > >
